|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 !****************************************************************************** 00007 MODULE powell 00008 USE kinds, ONLY: dp 00009 USE mathconstants, ONLY: twopi 00010 USE timings, ONLY: timeset,& 00011 timestop 00012 #include "cp_common_uses.h" 00013 00014 IMPLICIT NONE 00015 00016 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'powell' 00017 00018 TYPE opt_state_type 00019 INTEGER :: state 00020 INTEGER :: nvar 00021 INTEGER :: iprint 00022 INTEGER :: unit 00023 INTEGER :: maxfun 00024 REAL(dp) :: rhobeg, rhoend 00025 REAL(dp), DIMENSION(:), POINTER :: w 00026 REAL(dp), DIMENSION(:), POINTER :: xopt 00027 ! local variables 00028 INTEGER :: np, nh, nptm, nftest, idz, itest, nf, nfm, nfmm, 00029 nfsav, knew, kopt, ksave, ktemp 00030 REAL(dp) :: rhosq, recip, reciq, fbeg, fopt, diffa, xoptsq, 00031 rho, delta, dsq, dnorm, ratio, temp, tempq, beta, 00032 dx, vquad, diff, diffc, diffb, fsave, detrat, hdiag, 00033 distsq, gisq, gqsq, f, bstep, alpha, dstep 00034 END TYPE opt_state_type 00035 00036 PRIVATE 00037 PUBLIC :: powell_optimize, opt_state_type 00038 00039 !****************************************************************************** 00040 00041 CONTAINS 00042 00043 !****************************************************************************** 00044 SUBROUTINE powell_optimize (n,x,optstate) 00045 INTEGER :: n 00046 REAL(dp), DIMENSION(*) :: x 00047 TYPE(opt_state_type) :: optstate 00048 00049 CHARACTER(len=*), PARAMETER :: routineN = 'powell_optimize', 00050 routineP = moduleN//':'//routineN 00051 00052 INTEGER :: handle, npt 00053 00054 CALL timeset(routineN,handle) 00055 00056 SELECT CASE (optstate%state) 00057 CASE (0) 00058 npt=2*n+1 00059 ALLOCATE(optstate%w((npt+13)*(npt+n)+3*n*(n+3)/2)) 00060 ALLOCATE(optstate%xopt(n)) 00061 ! Initialize w 00062 optstate%w = 0.0_dp 00063 optstate%state = 1 00064 CALL newuoa (n, x, optstate) 00065 CASE (1,2) 00066 CALL newuoa (n, x, optstate) 00067 CASE (3) 00068 WRITE(optstate%unit,*) "POWELL| Exceeding maximum number of steps" 00069 optstate%state = -1 00070 CASE (4) 00071 WRITE(optstate%unit,*) "POWELL| Error in trust region" 00072 optstate%state = -1 00073 CASE (5) 00074 WRITE(optstate%unit,*) "POWELL| N out of range" 00075 optstate%state = -1 00076 CASE (6,7) 00077 optstate%state = -1 00078 CASE (8) 00079 x(1:n) = optstate%xopt(1:n) 00080 DEALLOCATE(optstate%w) 00081 DEALLOCATE(optstate%xopt) 00082 optstate%state = -1 00083 CASE DEFAULT 00084 STOP 00085 END SELECT 00086 00087 CALL timestop(handle) 00088 00089 END SUBROUTINE powell_optimize 00090 !****************************************************************************** 00091 SUBROUTINE newuoa (n,x,optstate) 00092 00093 INTEGER :: n 00094 REAL(dp), DIMENSION(*) :: x 00095 TYPE(opt_state_type) :: optstate 00096 00097 INTEGER :: ibmat, id, ifv, igq, ihq, 00098 ipq, ivl, iw, ixb, ixn, ixo, 00099 ixp, izmat, maxfun, ndim, np, 00100 npt, nptm 00101 REAL(dp) :: rhobeg, rhoend 00102 00103 maxfun = optstate%maxfun 00104 rhobeg = optstate%rhobeg 00105 rhoend = optstate%rhoend 00106 00107 ! 00108 ! This subroutine seeks the least value of a function of many variab 00109 ! by a trust region method that forms quadratic models by interpolat 00110 ! There can be some freedom in the interpolation conditions, which i 00111 ! taken up by minimizing the Frobenius norm of the change to the sec 00112 ! derivative of the quadratic model, beginning with a zero matrix. T 00113 ! arguments of the subroutine are as follows. 00114 ! 00115 ! N must be set to the number of variables and must be at least two. 00116 ! NPT is the number of interpolation conditions. Its value must be i 00117 ! interval [N+2,(N+1)(N+2)/2]. 00118 ! Initial values of the variables must be set in X(1),X(2),...,X(N). 00119 ! will be changed to the values that give the least calculated F. 00120 ! RHOBEG and RHOEND must be set to the initial and final values of a 00121 ! region radius, so both must be positive with RHOEND<=RHOBEG. Typ 00122 ! RHOBEG should be about one tenth of the greatest expected change 00123 ! variable, and RHOEND should indicate the accuracy that is requir 00124 ! the final values of the variables. 00125 ! The value of IPRINT should be set to 0, 1, 2 or 3, which controls 00126 ! amount of printing. Specifically, there is no output if IPRINT=0 00127 ! there is output only at the return if IPRINT=1. Otherwise, each 00128 ! value of RHO is printed, with the best vector of variables so fa 00129 ! the corresponding value of the objective function. Further, each 00130 ! value of F with its variables are output if IPRINT=3. 00131 ! MAXFUN must be set to an upper bound on the number of calls of CAL 00132 ! The array W will be used for working space. Its length must be at 00133 ! (NPT+13)*(NPT+N)+3*N*(N+3)/2. 00134 ! 00135 ! SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must se 00136 ! the value of the objective function for the variables X(1),X(2),.. 00137 ! 00138 ! Partition the working space array, so that different parts of it c 00139 ! treated separately by the subroutine that performs the main calcul 00140 ! 00141 np=n+1 00142 npt=2*n+1 00143 nptm=npt-np 00144 IF (npt < n+2 .OR. npt > ((n+2)*np)/2) THEN 00145 optstate%state = 5 00146 GO TO 20 00147 END IF 00148 ndim=npt+n 00149 ixb=1 00150 ixo=ixb+n 00151 ixn=ixo+n 00152 ixp=ixn+n 00153 ifv=ixp+n*npt 00154 igq=ifv+npt 00155 ihq=igq+n 00156 ipq=ihq+(n*np)/2 00157 ibmat=ipq+npt 00158 izmat=ibmat+ndim*n 00159 id=izmat+npt*nptm 00160 ivl=id+n 00161 iw=ivl+ndim 00162 ! 00163 ! The above settings provide a partition of W for subroutine NEWUOB. 00164 ! The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements 00165 ! W plus the space that is needed by the last array of NEWUOB. 00166 ! 00167 CALL newuob (n,npt,x,rhobeg,rhoend,maxfun,optstate%w(ixb:),optstate%w(ixo:),& 00168 optstate%w(ixn:),optstate%w(ixp:),optstate%w(ifv:),optstate%w(igq:),optstate%w(ihq:),& 00169 optstate%w(ipq:),optstate%w(ibmat:),optstate%w(izmat:),ndim,optstate%w(id:),& 00170 optstate%w(ivl:),optstate%w(iw:),optstate) 00171 00172 optstate%xopt(1:n) = optstate%w(ixb:ixb+n-1) + optstate%w(ixo:ixo+n-1) 00173 00174 20 CONTINUE 00175 00176 END SUBROUTINE newuoa 00177 00178 !****************************************************************************** 00179 SUBROUTINE newuob (n,npt,x,rhobeg,rhoend,maxfun,xbase,& 00180 xopt,xnew,xpt,fval,gq,hq,pq,bmat,zmat,ndim,d,vlag,w,opt) 00181 00182 INTEGER, INTENT(inout) :: n, npt 00183 REAL(dp), DIMENSION(1:n), INTENT(inout) :: x 00184 REAL(dp), INTENT(inout) :: rhobeg, rhoend 00185 INTEGER, INTENT(inout) :: maxfun 00186 REAL(dp), DIMENSION(*), INTENT(inout) :: xbase, xopt, xnew 00187 REAL(dp), DIMENSION(npt, *), 00188 INTENT(inout) :: xpt 00189 REAL(dp), DIMENSION(*), INTENT(inout) :: fval, gq, hq, pq 00190 INTEGER, INTENT(inout) :: ndim 00191 REAL(dp), DIMENSION(npt, *), 00192 INTENT(inout) :: zmat 00193 REAL(dp), DIMENSION(ndim, *), 00194 INTENT(inout) :: bmat 00195 REAL(dp), DIMENSION(*), INTENT(inout) :: d, vlag, w 00196 TYPE(opt_state_type) :: opt 00197 00198 INTEGER :: i, idz, ih, ip, ipt, itemp, 00199 itest, j, jp, jpt, k, knew, 00200 kopt, ksave, ktemp, nf, nfm, 00201 nfmm, nfsav, nftest, nh, np, 00202 nptm 00203 REAL(dp) :: alpha, beta, bstep, bsum, crvmin, delta, detrat, diff, diffa, 00204 diffb, diffc, distsq, dnorm, dsq, dstep, dx, f, fbeg, fopt, fsave, 00205 gisq, gqsq, half, hdiag, one, ratio, recip, reciq, rho, rhosq, sum, 00206 suma, sumb, sumz, temp, tempq, tenth, vquad, xipt, xjpt, xoptsq, zero 00207 00208 ! 00209 ! The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are ide 00210 ! to the corresponding arguments in SUBROUTINE NEWUOA. 00211 ! XBASE will hold a shift of origin that should reduce the contribut 00212 ! from rounding errors to values of the model and Lagrange functio 00213 ! XOPT will be set to the displacement from XBASE of the vector of 00214 ! variables that provides the least calculated F so far. 00215 ! XNEW will be set to the displacement from XBASE of the vector of 00216 ! variables for the current calculation of F. 00217 ! XPT will contain the interpolation point coordinates relative to X 00218 ! FVAL will hold the values of F at the interpolation points. 00219 ! GQ will hold the gradient of the quadratic model at XBASE. 00220 ! HQ will hold the explicit second derivatives of the quadratic mode 00221 ! PQ will contain the parameters of the implicit second derivatives 00222 ! the quadratic model. 00223 ! BMAT will hold the last N columns of H. 00224 ! ZMAT will hold the factorization of the leading NPT by NPT submatr 00225 ! H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, wh 00226 ! the elements of DZ are plus or minus one, as specified by IDZ. 00227 ! NDIM is the first dimension of BMAT and has the value NPT+N. 00228 ! D is reserved for trial steps from XOPT. 00229 ! VLAG will contain the values of the Lagrange functions at a new po 00230 ! They are part of a product that requires VLAG to be of length ND 00231 ! The array W will be used for working space. Its length must be at 00232 ! 10*NDIM = 10*(NPT+N). 00233 00234 IF ( opt%state == 1 ) THEN 00235 ! initialize all variable that will be stored 00236 np = 0 00237 nh = 0 00238 nptm = 0 00239 nftest = 0 00240 idz = 0 00241 itest = 0 00242 nf = 0 00243 nfm = 0 00244 nfmm = 0 00245 nfsav = 0 00246 knew = 0 00247 kopt = 0 00248 ksave = 0 00249 ktemp = 0 00250 rhosq = 0._dp 00251 recip = 0._dp 00252 reciq = 0._dp 00253 fbeg = 0._dp 00254 fopt = 0._dp 00255 diffa = 0._dp 00256 xoptsq = 0._dp 00257 rho = 0._dp 00258 delta = 0._dp 00259 dsq = 0._dp 00260 dnorm = 0._dp 00261 ratio = 0._dp 00262 temp = 0._dp 00263 tempq = 0._dp 00264 beta = 0._dp 00265 dx = 0._dp 00266 vquad = 0._dp 00267 diff = 0._dp 00268 diffc = 0._dp 00269 diffb = 0._dp 00270 fsave = 0._dp 00271 detrat = 0._dp 00272 hdiag = 0._dp 00273 distsq = 0._dp 00274 gisq = 0._dp 00275 gqsq = 0._dp 00276 f = 0._dp 00277 bstep = 0._dp 00278 alpha = 0._dp 00279 dstep = 0._dp 00280 ! 00281 END IF 00282 00283 ipt = 0 00284 jpt = 0 00285 xipt = 0._dp 00286 xjpt = 0._dp 00287 00288 half=0.5_dp 00289 one=1.0_dp 00290 tenth=0.1_dp 00291 zero=0.0_dp 00292 np=n+1 00293 nh=(n*np)/2 00294 nptm=npt-np 00295 nftest=MAX(maxfun,1) 00296 00297 IF ( opt%state == 2 ) GOTO 1000 00298 ! 00299 ! Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero. 00300 ! 00301 DO j=1,n 00302 xbase(j)=x(j) 00303 DO k=1,npt 00304 xpt(k,j)=zero 00305 END DO 00306 DO i=1,ndim 00307 bmat(i,j)=zero 00308 END DO 00309 END DO 00310 DO ih=1,nh 00311 hq(ih)=zero 00312 END DO 00313 DO k=1,npt 00314 pq(k)=zero 00315 DO j=1,nptm 00316 zmat(k,j)=zero 00317 END DO 00318 END DO 00319 ! 00320 ! Begin the initialization procedure. NF becomes one more than the n 00321 ! of function values so far. The coordinates of the displacement of 00322 ! next initial interpolation point from XBASE are set in XPT(NF,.). 00323 ! 00324 rhosq=rhobeg*rhobeg 00325 recip=one/rhosq 00326 reciq=SQRT(half)/rhosq 00327 nf=0 00328 50 nfm=nf 00329 nfmm=nf-n 00330 nf=nf+1 00331 IF (nfm <= 2*n) THEN 00332 IF (nfm >= 1 .AND. nfm <= N) THEN 00333 xpt(nf,nfm)=rhobeg 00334 ELSE IF (nfm > n) THEN 00335 xpt(nf,nfmm)=-rhobeg 00336 END IF 00337 ELSE 00338 itemp=(nfmm-1)/n 00339 jpt=nfm-itemp*n-n 00340 ipt=jpt+itemp 00341 IF (ipt > n) THEN 00342 itemp=jpt 00343 jpt=ipt-n 00344 ipt=itemp 00345 END IF 00346 xipt=rhobeg 00347 IF (fval(ipt+np) < fval(ipt+1)) xipt=-xipt 00348 XJPT=RHOBEG 00349 IF (fval(jpt+np) < fval(jpt+1)) xjpt=-xjpt 00350 xpt(nf,ipt)=xipt 00351 xpt(nf,jpt)=xjpt 00352 END IF 00353 ! 00354 ! Calculate the next value of F, label 70 being reached immediately 00355 ! after this calculation. The least function value so far and its in 00356 ! are required. 00357 ! 00358 DO j=1,n 00359 x(j)=xpt(nf,j)+xbase(j) 00360 END DO 00361 GOTO 310 00362 70 fval(nf)=f 00363 IF (nf == 1) THEN 00364 fbeg=f 00365 fopt=f 00366 kopt=1 00367 ELSE IF (f < fopt) THEN 00368 fopt=f 00369 kopt=nf 00370 END IF 00371 ! 00372 ! Set the nonzero initial elements of BMAT and the quadratic model i 00373 ! the cases when NF is at most 2*N+1. 00374 ! 00375 IF (NFM <= 2*N) THEN 00376 IF (nfm >= 1 .AND. nfm <= n) THEN 00377 gq(nfm)=(f-fbeg)/rhobeg 00378 IF (npt < nf+n) THEN 00379 bmat(1,nfm)=-one/rhobeg 00380 bmat(nf,nfm)=one/rhobeg 00381 bmat(npt+nfm,nfm)=-half*rhosq 00382 END IF 00383 ELSE IF (nfm > n) THEN 00384 bmat(nf-n,nfmm)=half/rhobeg 00385 bmat(nf,nfmm)=-half/rhobeg 00386 zmat(1,nfmm)=-reciq-reciq 00387 zmat(nf-n,nfmm)=reciq 00388 zmat(nf,nfmm)=reciq 00389 ih=(nfmm*(nfmm+1))/2 00390 temp=(fbeg-f)/rhobeg 00391 hq(ih)=(gq(nfmm)-temp)/rhobeg 00392 gq(nfmm)=half*(gq(nfmm)+temp) 00393 END IF 00394 ! 00395 ! Set the off-diagonal second derivatives of the Lagrange functions 00396 ! the initial quadratic model. 00397 ! 00398 ELSE 00399 ih=(ipt*(ipt-1))/2+jpt 00400 IF (xipt < zero) ipt=ipt+n 00401 IF (xjpt < zero) jpt=jpt+n 00402 zmat(1,nfmm)=recip 00403 zmat(nf,nfmm)=recip 00404 zmat(ipt+1,nfmm)=-recip 00405 zmat(jpt+1,nfmm)=-recip 00406 hq(ih)=(fbeg-fval(ipt+1)-fval(jpt+1)+f)/(xipt*xjpt) 00407 END IF 00408 IF (nf < npt) GOTO 50 00409 ! 00410 ! Begin the iterative procedure, because the initial model is comple 00411 ! 00412 rho=rhobeg 00413 delta=rho 00414 idz=1 00415 diffa=zero 00416 diffb=zero 00417 itest=0 00418 xoptsq=zero 00419 DO i=1,n 00420 xopt(i)=xpt(kopt,i) 00421 xoptsq=xoptsq+xopt(i)**2 00422 END DO 00423 90 nfsav=nf 00424 ! 00425 ! Generate the next trust region step and test its length. Set KNEW 00426 ! to -1 if the purpose of the next F will be to improve the model. 00427 ! 00428 100 knew=0 00429 CALL trsapp (n,npt,xopt,xpt,gq,hq,pq,delta,d,w,w(np),w(np+n),w(np+2*n),crvmin) 00430 dsq=zero 00431 DO i=1,n 00432 dsq=dsq+d(i)**2 00433 END DO 00434 dnorm=MIN(delta,SQRT(dsq)) 00435 IF (dnorm < half*rho) THEN 00436 knew=-1 00437 delta=tenth*delta 00438 ratio=-1.0_dp 00439 IF (delta <= 1.5_dp*rho) delta=rho 00440 IF (nf <= nfsav+2) GOTO 460 00441 temp=0.125_dp*crvmin*rho*rho 00442 IF (temp <= MAX(diffa,diffb,diffc)) GOTO 460 00443 GOTO 490 00444 END IF 00445 ! 00446 ! Shift XBASE if XOPT may be too far from XBASE. First make the chan 00447 ! to BMAT that do not depend on ZMAT. 00448 ! 00449 120 IF (dsq <= 1.0e-3_dp*xoptsq) THEN 00450 tempq=0.25_dp*xoptsq 00451 DO k=1,npt 00452 sum=zero 00453 DO i=1,n 00454 sum=sum+xpt(k,i)*xopt(i) 00455 END DO 00456 temp=pq(k)*sum 00457 sum=sum-half*xoptsq 00458 w(npt+k)=sum 00459 DO i=1,n 00460 gq(i)=gq(i)+temp*xpt(k,i) 00461 xpt(k,i)=xpt(k,i)-half*xopt(i) 00462 vlag(i)=bmat(k,i) 00463 w(i)=sum*xpt(k,i)+tempq*xopt(i) 00464 ip=npt+i 00465 DO j=1,i 00466 bmat(ip,j)=bmat(ip,j)+vlag(i)*w(j)+w(i)*vlag(j) 00467 END DO 00468 END DO 00469 END DO 00470 ! 00471 ! Then the revisions of BMAT that depend on ZMAT are calculated. 00472 ! 00473 DO k=1,nptm 00474 sumz=zero 00475 DO i=1,npt 00476 sumz=sumz+zmat(i,k) 00477 w(i)=w(npt+i)*zmat(i,k) 00478 END DO 00479 DO j=1,n 00480 sum=tempq*sumz*xopt(j) 00481 DO i=1,npt 00482 sum=sum+w(i)*xpt(i,j) 00483 vlag(j)=sum 00484 IF (k < idz) sum=-sum 00485 END DO 00486 DO i=1,npt 00487 bmat(i,j)=bmat(i,j)+sum*zmat(i,k) 00488 END DO 00489 END DO 00490 DO i=1,n 00491 ip=i+npt 00492 temp=vlag(i) 00493 IF (k < idz) temp=-temp 00494 DO j=1,i 00495 bmat(ip,j)=bmat(ip,j)+temp*vlag(j) 00496 END DO 00497 END DO 00498 END DO 00499 ! 00500 ! The following instructions complete the shift of XBASE, including 00501 ! the changes to the parameters of the quadratic model. 00502 ! 00503 ih=0 00504 DO j=1,n 00505 w(j)=zero 00506 DO k=1,npt 00507 w(j)=w(j)+pq(k)*xpt(k,j) 00508 xpt(k,j)=xpt(k,j)-half*xopt(j) 00509 END DO 00510 DO i=1,j 00511 ih=ih+1 00512 IF (i < j) gq(j)=gq(j)+hq(ih)*xopt(i) 00513 gq(i)=gq(i)+hq(ih)*xopt(j) 00514 hq(ih)=hq(ih)+w(i)*xopt(j)+xopt(i)*w(j) 00515 bmat(npt+i,j)=bmat(npt+j,i) 00516 END DO 00517 END DO 00518 DO j=1,n 00519 xbase(j)=xbase(j)+xopt(j) 00520 xopt(j)=zero 00521 END DO 00522 xoptsq=zero 00523 END IF 00524 ! 00525 ! Pick the model step if KNEW is positive. A different choice of D 00526 ! may be made later, if the choice of D by BIGLAG causes substantial 00527 ! cancellation in DENOM. 00528 ! 00529 IF (knew > 0) THEN 00530 CALL biglag (n,npt,xopt,xpt,bmat,zmat,idz,ndim,knew,dstep, & 00531 & d,alpha,vlag,vlag(npt+1),w,w(np),w(np+n)) 00532 END IF 00533 ! 00534 ! Calculate VLAG and BETA for the current choice of D. The first NPT 00535 ! components of W_check will be held in W. 00536 ! 00537 DO k=1,npt 00538 suma=zero 00539 sumb=zero 00540 sum=zero 00541 DO j=1,n 00542 suma=suma+xpt(k,j)*d(j) 00543 sumb=sumb+xpt(k,j)*xopt(j) 00544 sum=sum+bmat(k,j)*d(j) 00545 END DO 00546 w(k)=suma*(half*suma+sumb) 00547 vlag(k)=sum 00548 END DO 00549 beta=zero 00550 DO k=1,nptm 00551 sum=zero 00552 DO i=1,npt 00553 sum=sum+zmat(i,k)*w(i) 00554 END DO 00555 IF (k < idz) THEN 00556 beta=beta+sum*sum 00557 sum=-sum 00558 ELSE 00559 beta=beta-sum*sum 00560 END IF 00561 DO i=1,npt 00562 vlag(i)=vlag(i)+sum*zmat(i,k) 00563 END DO 00564 END DO 00565 bsum=zero 00566 dx=zero 00567 DO j=1,n 00568 sum=zero 00569 DO i=1,npt 00570 sum=sum+w(i)*bmat(i,j) 00571 END DO 00572 bsum=bsum+sum*d(j) 00573 jp=npt+j 00574 DO k=1,n 00575 sum=sum+bmat(jp,k)*d(k) 00576 END DO 00577 vlag(jp)=sum 00578 bsum=bsum+sum*d(j) 00579 dx=dx+d(j)*xopt(j) 00580 END DO 00581 beta=dx*dx+dsq*(xoptsq+dx+dx+half*dsq)+beta-bsum 00582 vlag(kopt)=vlag(kopt)+one 00583 ! 00584 ! If KNEW is positive and if the cancellation in DENOM is unacceptab 00585 ! then BIGDEN calculates an alternative model step, XNEW being used 00586 ! working space. 00587 ! 00588 IF (knew > 0) THEN 00589 temp=one+alpha*beta/vlag(knew)**2 00590 IF (ABS(temp) <= 0.8_dp) THEN 00591 CALL bigden (n,npt,xopt,xpt,bmat,zmat,idz,ndim,kopt, & 00592 & knew,d,w,vlag,beta,xnew,w(ndim+1),w(6*ndim+1)) 00593 END IF 00594 END IF 00595 ! 00596 ! Calculate the next value of the objective function. 00597 ! 00598 290 DO i=1,n 00599 xnew(i)=xopt(i)+d(i) 00600 x(i)=xbase(i)+xnew(i) 00601 END DO 00602 nf=nf+1 00603 310 IF (nf > nftest) THEN 00604 ! return to many steps 00605 nf=nf-1 00606 opt%state = 3 00607 CALL get_state 00608 GOTO 530 00609 END IF 00610 00611 CALL get_state 00612 00613 opt%state = 2 00614 00615 RETURN 00616 00617 1000 CONTINUE 00618 00619 CALL set_state 00620 00621 IF (nf <= npt) GOTO 70 00622 IF (knew == -1) THEN 00623 opt%state = 6 00624 CALL get_state 00625 GOTO 530 00626 END IF 00627 ! 00628 ! Use the quadratic model to predict the change in F due to the step 00629 ! and set DIFF to the error of this prediction. 00630 ! 00631 vquad=zero 00632 ih=0 00633 DO j=1,n 00634 vquad=vquad+d(j)*gq(j) 00635 DO i=1,j 00636 ih=ih+1 00637 temp=d(i)*xnew(j)+d(j)*xopt(i) 00638 IF (i == j) temp=half*temp 00639 vquad=vquad+temp*hq(ih) 00640 END DO 00641 END DO 00642 DO k=1,npt 00643 vquad=vquad+pq(k)*w(k) 00644 END DO 00645 diff=f-fopt-vquad 00646 diffc=diffb 00647 diffb=diffa 00648 diffa=ABS(diff) 00649 IF (dnorm > rho) nfsav=nf 00650 ! 00651 ! Update FOPT and XOPT if the new F is the least value of the object 00652 ! function so far. The branch when KNEW is positive occurs if D is n 00653 ! a trust region step. 00654 ! 00655 fsave=fopt 00656 IF (f < fopt) THEN 00657 fopt=f 00658 xoptsq=zero 00659 DO i=1,n 00660 xopt(i)=xnew(i) 00661 xoptsq=xoptsq+xopt(i)**2 00662 END DO 00663 END IF 00664 ksave=knew 00665 IF (knew > 0) GOTO 410 00666 ! 00667 ! Pick the next value of DELTA after a trust region step. 00668 ! 00669 IF (vquad >= zero) THEN 00670 ! Return because a trust region step has failed to reduce Q 00671 opt%state = 4 00672 CALL get_state 00673 GOTO 530 00674 END IF 00675 ratio=(f-fsave)/vquad 00676 IF (ratio <= tenth) THEN 00677 delta=half*dnorm 00678 ELSE IF (ratio <= 0.7_dp) THEN 00679 delta=MAX(half*delta,dnorm) 00680 ELSE 00681 delta=MAX(half*delta,dnorm+dnorm) 00682 END IF 00683 IF (delta <= 1.5_dp*rho) delta=rho 00684 ! 00685 ! Set KNEW to the index of the next interpolation point to be delete 00686 ! 00687 rhosq=MAX(tenth*delta,rho)**2 00688 ktemp=0 00689 detrat=zero 00690 IF (f >= fsave) THEN 00691 ktemp=kopt 00692 detrat=one 00693 END IF 00694 DO k=1,npt 00695 hdiag=zero 00696 DO j=1,nptm 00697 temp=one 00698 IF (j < idz) temp=-one 00699 hdiag=hdiag+temp*zmat(k,j)**2 00700 END DO 00701 temp=ABS(beta*hdiag+vlag(k)**2) 00702 distsq=zero 00703 DO j=1,n 00704 distsq=distsq+(xpt(k,j)-xopt(j))**2 00705 END DO 00706 IF (distsq > rhosq) temp=temp*(distsq/rhosq)**3 00707 IF (temp > detrat .AND. k /= ktemp) THEN 00708 detrat=temp 00709 knew=k 00710 END IF 00711 END DO 00712 IF (knew == 0) GOTO 460 00713 ! 00714 ! Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point 00715 ! can be moved. Begin the updating of the quadratic model, starting 00716 ! with the explicit second derivative term. 00717 ! 00718 410 CALL update (n,npt,bmat,zmat,idz,ndim,vlag,beta,knew,w) 00719 fval(knew)=f 00720 ih=0 00721 DO i=1,n 00722 temp=pq(knew)*xpt(knew,i) 00723 DO j=1,i 00724 ih=ih+1 00725 hq(ih)=hq(ih)+temp*xpt(knew,j) 00726 END DO 00727 END DO 00728 pq(knew)=zero 00729 ! 00730 ! Update the other second derivative parameters, and then the gradie 00731 ! vector of the model. Also include the new interpolation point. 00732 ! 00733 DO j=1,nptm 00734 temp=diff*zmat(knew,j) 00735 IF (j < idz) temp=-temp 00736 DO k=1,npt 00737 pq(k)=pq(k)+temp*zmat(k,j) 00738 END DO 00739 END DO 00740 gqsq=zero 00741 DO i=1,n 00742 gq(i)=gq(i)+diff*bmat(knew,i) 00743 gqsq=gqsq+gq(i)**2 00744 xpt(knew,i)=xnew(i) 00745 END DO 00746 ! 00747 ! If a trust region step makes a small change to the objective funct 00748 ! then calculate the gradient of the least Frobenius norm interpolan 00749 ! XBASE, and store it in W, using VLAG for a vector of right hand si 00750 ! 00751 IF (ksave == 0 .AND. delta == rho) THEN 00752 IF (ABS(ratio) > 1.0e-2_dp) THEN 00753 itest=0 00754 ELSE 00755 DO k=1,npt 00756 vlag(k)=fval(k)-fval(kopt) 00757 END DO 00758 gisq=zero 00759 DO i=1,n 00760 sum=zero 00761 DO k=1,npt 00762 sum=sum+bmat(k,i)*vlag(k) 00763 END DO 00764 gisq=gisq+sum*sum 00765 w(i)=sum 00766 END DO 00767 ! 00768 ! Test whether to replace the new quadratic model by the least Frobe 00769 ! norm interpolant, making the replacement if the test is satisfied. 00770 ! 00771 itest=itest+1 00772 IF (gqsq < 1.0e2_dp*gisq) itest=0 00773 IF (itest >= 3) THEN 00774 DO i=1,n 00775 gq(i)=w(i) 00776 END DO 00777 DO ih=1,nh 00778 hq(ih)=zero 00779 END DO 00780 DO j=1,nptm 00781 w(j)=zero 00782 DO k=1,npt 00783 w(j)=w(j)+vlag(k)*zmat(k,j) 00784 END DO 00785 IF (j < idz) w(j)=-w(j) 00786 END DO 00787 DO k=1,npt 00788 pq(k)=zero 00789 DO j=1,nptm 00790 pq(k)=pq(k)+zmat(k,j)*w(j) 00791 END DO 00792 END DO 00793 itest=0 00794 END IF 00795 END IF 00796 END IF 00797 IF (f < fsave) kopt=knew 00798 ! 00799 ! If a trust region step has provided a sufficient decrease in F, th 00800 ! branch for another trust region calculation. The case KSAVE>0 occu 00801 ! when the new function value was calculated by a model step. 00802 ! 00803 IF (f <= fsave+tenth*vquad) GOTO 100 00804 IF (ksave > 0) GOTO 100 00805 ! 00806 ! Alternatively, find out if the interpolation points are close enou 00807 ! to the best point so far. 00808 ! 00809 knew=0 00810 460 distsq=4.0_dp*delta*delta 00811 DO k=1,npt 00812 sum=zero 00813 DO j=1,n 00814 sum=sum+(xpt(k,j)-xopt(j))**2 00815 END DO 00816 IF (sum > distsq) THEN 00817 knew=k 00818 distsq=sum 00819 END IF 00820 END DO 00821 ! 00822 ! If KNEW is positive, then set DSTEP, and branch back for the next 00823 ! iteration, which will generate a "model step". 00824 ! 00825 IF (knew > 0) THEN 00826 dstep=MAX(MIN(tenth*SQRT(distsq),half*delta),rho) 00827 dsq=dstep*dstep 00828 GOTO 120 00829 END IF 00830 IF (ratio > zero) GOTO 100 00831 IF (MAX(delta,dnorm) > rho) GOTO 100 00832 ! 00833 ! The calculations with the current value of RHO are complete. Pick 00834 ! next values of RHO and DELTA. 00835 ! 00836 490 IF (rho > rhoend) THEN 00837 delta=half*rho 00838 ratio=rho/rhoend 00839 IF (ratio <= 16.0_dp) THEN 00840 rho=rhoend 00841 ELSE IF (ratio <= 250.0_dp) THEN 00842 rho=SQRT(ratio)*rhoend 00843 ELSE 00844 rho=tenth*rho 00845 END IF 00846 delta=MAX(delta,rho) 00847 GOTO 90 00848 END IF 00849 ! 00850 ! Return from the calculation, after another Newton-Raphson step, if 00851 ! it is too short to have been tried before. 00852 ! 00853 IF (knew == -1) GOTO 290 00854 opt%state = 7 00855 CALL get_state 00856 530 IF (fopt <= f) THEN 00857 DO i=1,n 00858 x(i)=xbase(i)+xopt(i) 00859 END DO 00860 f=fopt 00861 END IF 00862 00863 CALL get_state 00864 00865 !****************************************************************************** 00866 CONTAINS 00867 !****************************************************************************** 00868 SUBROUTINE get_state 00869 opt%np = np 00870 opt%nh = nh 00871 opt%nptm = nptm 00872 opt%nftest = nftest 00873 opt%idz = idz 00874 opt%itest = itest 00875 opt%nf = nf 00876 opt%nfm = nfm 00877 opt%nfmm = nfmm 00878 opt%nfsav = nfsav 00879 opt%knew = knew 00880 opt%kopt = kopt 00881 opt%ksave = ksave 00882 opt%ktemp = ktemp 00883 opt%rhosq = rhosq 00884 opt%recip = recip 00885 opt%reciq = reciq 00886 opt%fbeg = fbeg 00887 opt%fopt = fopt 00888 opt%diffa = diffa 00889 opt%xoptsq = xoptsq 00890 opt%rho = rho 00891 opt%delta = delta 00892 opt%dsq = dsq 00893 opt%dnorm = dnorm 00894 opt%ratio = ratio 00895 opt%temp = temp 00896 opt%tempq = tempq 00897 opt%beta = beta 00898 opt%dx = dx 00899 opt%vquad = vquad 00900 opt%diff = diff 00901 opt%diffc = diffc 00902 opt%diffb = diffb 00903 opt%fsave = fsave 00904 opt%detrat = detrat 00905 opt%hdiag = hdiag 00906 opt%distsq = distsq 00907 opt%gisq = gisq 00908 opt%gqsq = gqsq 00909 opt%f = f 00910 opt%bstep = bstep 00911 opt%alpha = alpha 00912 opt%dstep = dstep 00913 END SUBROUTINE get_state 00914 00915 !****************************************************************************** 00916 SUBROUTINE set_state 00917 np = opt%np 00918 nh = opt%nh 00919 nptm = opt%nptm 00920 nftest = opt%nftest 00921 idz = opt%idz 00922 itest = opt%itest 00923 nf = opt%nf 00924 nfm = opt%nfm 00925 nfmm = opt%nfmm 00926 nfsav = opt%nfsav 00927 knew = opt%knew 00928 kopt = opt%kopt 00929 ksave = opt%ksave 00930 ktemp = opt%ktemp 00931 rhosq = opt%rhosq 00932 recip = opt%recip 00933 reciq = opt%reciq 00934 fbeg = opt%fbeg 00935 fopt = opt%fopt 00936 diffa = opt%diffa 00937 xoptsq = opt%xoptsq 00938 rho = opt%rho 00939 delta = opt%delta 00940 dsq = opt%dsq 00941 dnorm = opt%dnorm 00942 ratio = opt%ratio 00943 temp = opt%temp 00944 tempq = opt%tempq 00945 beta = opt%beta 00946 dx = opt%dx 00947 vquad = opt%vquad 00948 diff = opt%diff 00949 diffc = opt%diffc 00950 diffb = opt%diffb 00951 fsave = opt%fsave 00952 detrat = opt%detrat 00953 hdiag = opt%hdiag 00954 distsq = opt%distsq 00955 gisq = opt%gisq 00956 gqsq = opt%gqsq 00957 f = opt%f 00958 bstep = opt%bstep 00959 alpha = opt%alpha 00960 dstep = opt%dstep 00961 END SUBROUTINE set_state 00962 00963 END SUBROUTINE newuob 00964 00965 !****************************************************************************** 00966 00967 SUBROUTINE bigden (n,npt,xopt,xpt,bmat,zmat,idz,ndim,kopt,& 00968 knew,d,w,vlag,beta,s,wvec,prod) 00969 00970 INTEGER, INTENT(inout) :: n, npt 00971 REAL(dp), DIMENSION(*), INTENT(inout) :: xopt 00972 REAL(dp), DIMENSION(npt, *), 00973 INTENT(inout) :: xpt 00974 INTEGER, INTENT(inout) :: ndim, idz 00975 REAL(dp), DIMENSION(npt, *), 00976 INTENT(inout) :: zmat 00977 REAL(dp), DIMENSION(ndim, *), 00978 INTENT(inout) :: bmat 00979 INTEGER, INTENT(inout) :: kopt, knew 00980 REAL(dp), DIMENSION(*), INTENT(inout) :: d, w, vlag 00981 REAL(dp), INTENT(inout) :: beta 00982 REAL(dp), DIMENSION(*), INTENT(inout) :: s 00983 REAL(dp), DIMENSION(ndim, *), 00984 INTENT(inout) :: wvec, prod 00985 00986 REAL(dp), PARAMETER :: half = 0.5_dp, one = 1._dp, 00987 quart = 0.25_dp, two = 2._dp, 00988 zero = 0._dp 00989 00990 INTEGER :: i, ip, isave, iterc, iu, j, 00991 jc, k, ksav, nptm, nw 00992 REAL(dp) :: alpha, angle, dd, denmax, denold, densav, diff, ds, dstemp, 00993 dtest, ss, ssden, sstemp, step, sum, sumold, tau, temp, tempa, tempb, 00994 tempc, xoptd, xopts, xoptsq 00995 REAL(dp), DIMENSION(9) :: den, denex, par 00996 00997 ! 00998 ! N is the number of variables. 00999 ! NPT is the number of interpolation equations. 01000 ! XOPT is the best interpolation point so far. 01001 ! XPT contains the coordinates of the current interpolation points. 01002 ! BMAT provides the last N columns of H. 01003 ! ZMAT and IDZ give a factorization of the first NPT by NPT submatri 01004 ! NDIM is the first dimension of BMAT and has the value NPT+N. 01005 ! KOPT is the index of the optimal interpolation point. 01006 ! KNEW is the index of the interpolation point that is going to be m 01007 ! D will be set to the step from XOPT to the new point, and on entry 01008 ! should be the D that was calculated by the last call of BIGLAG. 01009 ! length of the initial D provides a trust region bound on the fin 01010 ! W will be set to Wcheck for the final choice of D. 01011 ! VLAG will be set to Theta*Wcheck+e_b for the final choice of D. 01012 ! BETA will be set to the value that will occur in the updating form 01013 ! when the KNEW-th interpolation point is moved to its new positio 01014 ! S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be us 01015 ! for working space. 01016 ! 01017 ! D is calculated in a way that should provide a denominator with a 01018 ! modulus in the updating formula when the KNEW-th interpolation poi 01019 ! shifted to the new position XOPT+D. 01020 ! 01021 01022 nptm=npt-n-1 01023 ! 01024 ! Store the first NPT elements of the KNEW-th column of H in W(N+1) 01025 ! to W(N+NPT). 01026 ! 01027 DO k=1,npt 01028 w(n+k)=zero 01029 END DO 01030 DO j=1,nptm 01031 temp=zmat(knew,j) 01032 IF (j < idz) temp=-temp 01033 DO k=1,npt 01034 w(n+k)=w(n+k)+temp*zmat(k,j) 01035 END DO 01036 END DO 01037 alpha=w(n+knew) 01038 ! 01039 ! The initial search direction D is taken from the last call of BIGL 01040 ! and the initial S is set below, usually to the direction from X_OP 01041 ! to X_KNEW, but a different direction to an interpolation point may 01042 ! be chosen, in order to prevent S from being nearly parallel to D. 01043 ! 01044 dd=zero 01045 ds=zero 01046 ss=zero 01047 xoptsq=zero 01048 DO i=1,n 01049 dd=dd+d(i)**2 01050 s(i)=xpt(knew,i)-xopt(i) 01051 ds=ds+d(i)*s(i) 01052 ss=ss+s(i)**2 01053 xoptsq=xoptsq+xopt(i)**2 01054 END DO 01055 IF (ds*ds > 0.99_dp*dd*ss) THEN 01056 ksav=knew 01057 dtest=ds*ds/ss 01058 DO k=1,npt 01059 IF (k /= kopt) THEN 01060 dstemp=zero 01061 sstemp=zero 01062 DO i=1,n 01063 diff=xpt(k,i)-xopt(i) 01064 dstemp=dstemp+d(i)*diff 01065 sstemp=sstemp+diff*diff 01066 END DO 01067 IF (dstemp*dstemp/sstemp < dtest) THEN 01068 ksav=k 01069 dtest=dstemp*dstemp/sstemp 01070 ds=dstemp 01071 ss=sstemp 01072 END IF 01073 END IF 01074 END DO 01075 DO i=1,n 01076 s(i)=xpt(ksav,i)-xopt(i) 01077 END DO 01078 END IF 01079 ssden=dd*ss-ds*ds 01080 iterc=0 01081 densav=zero 01082 ! 01083 ! Begin the iteration by overwriting S with a vector that has the 01084 ! required length and direction. 01085 ! 01086 mainloop : DO 01087 iterc=iterc+1 01088 temp=one/SQRT(ssden) 01089 xoptd=zero 01090 xopts=zero 01091 DO i=1,n 01092 s(i)=temp*(dd*s(i)-ds*d(i)) 01093 xoptd=xoptd+xopt(i)*d(i) 01094 xopts=xopts+xopt(i)*s(i) 01095 END DO 01096 ! 01097 ! Set the coefficients of the first two terms of BETA. 01098 ! 01099 tempa=half*xoptd*xoptd 01100 tempb=half*xopts*xopts 01101 den(1)=dd*(xoptsq+half*dd)+tempa+tempb 01102 den(2)=two*xoptd*dd 01103 den(3)=two*xopts*dd 01104 den(4)=tempa-tempb 01105 den(5)=xoptd*xopts 01106 DO i=6,9 01107 den(i)=zero 01108 END DO 01109 ! 01110 ! Put the coefficients of Wcheck in WVEC. 01111 ! 01112 DO k=1,npt 01113 tempa=zero 01114 tempb=zero 01115 tempc=zero 01116 DO i=1,n 01117 tempa=tempa+xpt(k,i)*d(i) 01118 tempb=tempb+xpt(k,i)*s(i) 01119 tempc=tempc+xpt(k,i)*xopt(i) 01120 END DO 01121 wvec(k,1)=quart*(tempa*tempa+tempb*tempb) 01122 wvec(k,2)=tempa*tempc 01123 wvec(k,3)=tempb*tempc 01124 wvec(k,4)=quart*(tempa*tempa-tempb*tempb) 01125 wvec(k,5)=half*tempa*tempb 01126 END DO 01127 DO i=1,n 01128 ip=i+npt 01129 wvec(ip,1)=zero 01130 wvec(ip,2)=d(i) 01131 wvec(ip,3)=s(i) 01132 wvec(ip,4)=zero 01133 wvec(ip,5)=zero 01134 END DO 01135 ! 01136 ! Put the coefficents of THETA*Wcheck in PROD. 01137 ! 01138 DO jc=1,5 01139 nw=npt 01140 IF (jc == 2 .OR. jc == 3) nw=ndim 01141 DO k=1,npt 01142 prod(k,jc)=zero 01143 END DO 01144 DO j=1,nptm 01145 sum=zero 01146 DO k=1,npt 01147 sum=sum+zmat(k,j)*wvec(k,jc) 01148 END DO 01149 IF (j < idz) sum=-sum 01150 DO k=1,npt 01151 prod(k,jc)=prod(k,jc)+sum*zmat(k,j) 01152 END DO 01153 END DO 01154 IF (nw == ndim) THEN 01155 DO k=1,npt 01156 sum=zero 01157 DO j=1,n 01158 sum=sum+bmat(k,j)*wvec(npt+j,jc) 01159 END DO 01160 prod(k,jc)=prod(k,jc)+sum 01161 END DO 01162 END IF 01163 DO j=1,n 01164 sum=zero 01165 DO i=1,nw 01166 sum=sum+bmat(i,j)*wvec(i,jc) 01167 END DO 01168 prod(npt+j,jc)=sum 01169 END DO 01170 END DO 01171 ! 01172 ! Include in DEN the part of BETA that depends on THETA. 01173 ! 01174 DO k=1,ndim 01175 sum=zero 01176 DO I=1,5 01177 par(i)=half*prod(k,i)*wvec(k,i) 01178 sum=sum+par(i) 01179 END DO 01180 den(1)=den(1)-par(1)-sum 01181 tempa=prod(k,1)*wvec(k,2)+prod(k,2)*wvec(k,1) 01182 tempb=prod(k,2)*wvec(k,4)+prod(k,4)*wvec(k,2) 01183 tempc=prod(k,3)*wvec(k,5)+prod(k,5)*wvec(k,3) 01184 den(2)=den(2)-tempa-half*(tempb+tempc) 01185 den(6)=den(6)-half*(tempb-tempc) 01186 tempa=prod(k,1)*wvec(k,3)+prod(k,3)*wvec(k,1) 01187 tempb=prod(k,2)*wvec(k,5)+prod(k,5)*wvec(k,2) 01188 tempc=prod(k,3)*wvec(k,4)+prod(k,4)*wvec(k,3) 01189 den(3)=den(3)-tempa-half*(tempb-tempc) 01190 den(7)=den(7)-half*(tempb+tempc) 01191 tempa=prod(k,1)*wvec(k,4)+prod(k,4)*wvec(k,1) 01192 den(4)=den(4)-tempa-par(2)+par(3) 01193 tempa=prod(k,1)*wvec(k,5)+prod(k,5)*wvec(k,1) 01194 tempb=prod(k,2)*wvec(k,3)+prod(k,3)*wvec(k,2) 01195 den(5)=den(5)-tempa-half*tempb 01196 den(8)=den(8)-par(4)+par(5) 01197 tempa=prod(k,4)*wvec(k,5)+prod(k,5)*wvec(k,4) 01198 den(9)=den(9)-half*tempa 01199 END DO 01200 ! 01201 ! Extend DEN so that it holds all the coefficients of DENOM. 01202 ! 01203 sum=zero 01204 DO i=1,5 01205 par(i)=half*prod(knew,i)**2 01206 sum=sum+par(i) 01207 END DO 01208 denex(1)=alpha*den(1)+par(1)+sum 01209 tempa=two*prod(knew,1)*prod(knew,2) 01210 tempb=prod(knew,2)*prod(knew,4) 01211 tempc=prod(knew,3)*prod(knew,5) 01212 denex(2)=alpha*den(2)+tempa+tempb+tempc 01213 denex(6)=alpha*den(6)+tempb-tempc 01214 tempa=two*prod(knew,1)*prod(knew,3) 01215 tempb=prod(knew,2)*prod(knew,5) 01216 tempc=prod(knew,3)*prod(knew,4) 01217 denex(3)=alpha*den(3)+tempa+tempb-tempc 01218 denex(7)=alpha*den(7)+tempb+tempc 01219 tempa=two*prod(knew,1)*prod(knew,4) 01220 denex(4)=alpha*den(4)+tempa+par(2)-par(3) 01221 tempa=two*prod(knew,1)*prod(knew,5) 01222 denex(5)=alpha*den(5)+tempa+prod(knew,2)*prod(knew,3) 01223 denex(8)=alpha*den(8)+par(4)-par(5) 01224 denex(9)=alpha*den(9)+prod(knew,4)*prod(knew,5) 01225 ! 01226 ! Seek the value of the angle that maximizes the modulus of DENOM. 01227 ! 01228 sum=denex(1)+denex(2)+denex(4)+denex(6)+denex(8) 01229 denold=sum 01230 denmax=sum 01231 isave=0 01232 iu=49 01233 temp=twopi/REAL(IU+1,dp) 01234 par(1)=one 01235 DO i=1,iu 01236 angle=REAL(i,dp)*temp 01237 par(2)=COS(angle) 01238 par(3)=SIN(angle) 01239 DO j=4,8,2 01240 par(j)=par(2)*par(j-2)-par(3)*par(j-1) 01241 par(j+1)=par(2)*par(j-1)+par(3)*par(j-2) 01242 END DO 01243 sumold=sum 01244 sum=zero 01245 DO j=1,9 01246 sum=sum+denex(j)*par(j) 01247 END DO 01248 IF (ABS(sum) > ABS(denmax)) THEN 01249 denmax=sum 01250 isave=i 01251 tempa=sumold 01252 ELSE IF (i == isave+1) THEN 01253 tempb=sum 01254 END IF 01255 END DO 01256 IF (isave == 0) tempa=sum 01257 IF (isave == iu) tempb=denold 01258 step=zero 01259 IF (tempa /= tempb) THEN 01260 tempa=tempa-denmax 01261 tempb=tempb-denmax 01262 step=half*(tempa-tempb)/(tempa+tempb) 01263 END IF 01264 angle=temp*(REAL(isave,dp)+step) 01265 ! 01266 ! Calculate the new parameters of the denominator, the new VLAG vect 01267 ! and the new D. Then test for convergence. 01268 ! 01269 par(2)=COS(angle) 01270 par(3)=SIN(angle) 01271 DO j=4,8,2 01272 par(j)=par(2)*par(j-2)-par(3)*par(j-1) 01273 par(j+1)=par(2)*par(j-1)+par(3)*par(j-2) 01274 END DO 01275 beta=zero 01276 denmax=zero 01277 DO j=1,9 01278 beta=beta+den(j)*par(j) 01279 denmax=denmax+denex(j)*par(j) 01280 END DO 01281 DO k=1,ndim 01282 vlag(k)=zero 01283 DO j=1,5 01284 vlag(k)=vlag(k)+prod(k,j)*par(j) 01285 END DO 01286 END DO 01287 tau=vlag(knew) 01288 dd=zero 01289 tempa=zero 01290 tempb=zero 01291 DO i=1,n 01292 d(i)=par(2)*d(i)+par(3)*s(i) 01293 w(i)=xopt(i)+d(i) 01294 dd=dd+d(i)**2 01295 tempa=tempa+d(i)*w(i) 01296 tempb=tempb+w(i)*w(i) 01297 END DO 01298 IF (iterc >= n) EXIT mainloop 01299 IF (iterc >= 1) densav=MAX(densav,denold) 01300 IF (ABS(denmax) <= 1.1_dp*ABS(densav)) EXIT mainloop 01301 densav=denmax 01302 ! 01303 ! Set S to half the gradient of the denominator with respect to D. 01304 ! Then branch for the next iteration. 01305 ! 01306 DO i=1,n 01307 temp=tempa*xopt(i)+tempb*d(i)-vlag(npt+i) 01308 s(i)=tau*bmat(knew,i)+alpha*temp 01309 END DO 01310 DO k=1,npt 01311 sum=zero 01312 DO j=1,n 01313 sum=sum+xpt(k,j)*w(j) 01314 END DO 01315 temp=(tau*w(n+k)-alpha*vlag(k))*sum 01316 DO i=1,n 01317 s(i)=s(i)+temp*xpt(k,i) 01318 END DO 01319 END DO 01320 ss=zero 01321 ds=zero 01322 DO i=1,n 01323 ss=ss+s(i)**2 01324 ds=ds+d(i)*s(i) 01325 END DO 01326 ssden=dd*ss-ds*ds 01327 IF (ssden < 1.0e-8_dp*dd*ss) EXIT mainloop 01328 END DO mainloop 01329 ! 01330 ! Set the vector W before the RETURN from the subroutine. 01331 ! 01332 DO k=1,ndim 01333 w(k)=zero 01334 DO j=1,5 01335 w(k)=w(k)+wvec(k,j)*par(j) 01336 END DO 01337 END DO 01338 vlag(kopt)=vlag(kopt)+one 01339 01340 END SUBROUTINE bigden 01341 !****************************************************************************** 01342 01343 SUBROUTINE biglag (n,npt,xopt,xpt,bmat,zmat,idz,ndim,knew,& 01344 delta,d,alpha,hcol,gc,gd,s,w) 01345 INTEGER, INTENT(inout) :: n, npt 01346 REAL(dp), DIMENSION(*), INTENT(inout) :: xopt 01347 REAL(dp), DIMENSION(npt, *), 01348 INTENT(inout) :: xpt 01349 INTEGER, INTENT(inout) :: ndim, idz 01350 REAL(dp), DIMENSION(npt, *), 01351 INTENT(inout) :: zmat 01352 REAL(dp), DIMENSION(ndim, *), 01353 INTENT(inout) :: bmat 01354 INTEGER, INTENT(inout) :: knew 01355 REAL(dp), INTENT(inout) :: delta 01356 REAL(dp), DIMENSION(*), INTENT(inout) :: d 01357 REAL(dp), INTENT(inout) :: alpha 01358 REAL(dp), DIMENSION(*), INTENT(inout) :: hcol, gc, gd, s, w 01359 01360 REAL(dp), PARAMETER :: half = 0.5_dp, one = 1._dp, 01361 zero = 0._dp 01362 01363 INTEGER :: i, isave, iterc, iu, j, k, 01364 nptm 01365 REAL(dp) :: angle, cf1, cf2, cf3, cf4, cf5, cth, dd, delsq, denom, dhd, 01366 gg, scale, sp, ss, step, sth, sum, tau, taubeg, taumax, tauold, temp, 01367 tempa, tempb 01368 01369 ! 01370 ! 01371 ! N is the number of variables. 01372 ! NPT is the number of interpolation equations. 01373 ! XOPT is the best interpolation point so far. 01374 ! XPT contains the coordinates of the current interpolation points. 01375 ! BMAT provides the last N columns of H. 01376 ! ZMAT and IDZ give a factorization of the first NPT by NPT submatrix 01377 ! NDIM is the first dimension of BMAT and has the value NPT+N. 01378 ! KNEW is the index of the interpolation point that is going to be m 01379 ! DELTA is the current trust region bound. 01380 ! D will be set to the step from XOPT to the new point. 01381 ! ALPHA will be set to the KNEW-th diagonal element of the H matrix. 01382 ! HCOL, GC, GD, S and W will be used for working space. 01383 ! 01384 ! The step D is calculated in a way that attempts to maximize the mo 01385 ! of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFU 01386 ! the KNEW-th Lagrange function. 01387 ! 01388 01389 delsq=delta*delta 01390 nptm=npt-n-1 01391 ! 01392 ! Set the first NPT components of HCOL to the leading elements of th 01393 ! KNEW-th column of H. 01394 ! 01395 iterc=0 01396 DO k=1,npt 01397 hcol(k)=zero 01398 END DO 01399 DO j=1,nptm 01400 temp=zmat(knew,j) 01401 IF (j < idz) temp=-temp 01402 DO k=1,npt 01403 hcol(k)=hcol(k)+temp*zmat(k,j) 01404 END DO 01405 END DO 01406 alpha=hcol(knew) 01407 ! 01408 ! Set the unscaled initial direction D. Form the gradient of LFUNC a 01409 ! XOPT, and multiply D by the second derivative matrix of LFUNC. 01410 ! 01411 dd=zero 01412 DO i=1,n 01413 d(i)=xpt(knew,i)-xopt(i) 01414 gc(i)=bmat(knew,i) 01415 gd(i)=zero 01416 dd=dd+d(i)**2 01417 END DO 01418 DO k=1,npt 01419 temp=zero 01420 sum=zero 01421 DO j=1,n 01422 temp=temp+xpt(k,j)*xopt(j) 01423 sum=sum+xpt(k,j)*d(j) 01424 END DO 01425 temp=hcol(k)*temp 01426 sum=hcol(k)*sum 01427 DO i=1,n 01428 gc(i)=gc(i)+temp*xpt(k,i) 01429 gd(i)=gd(i)+sum*xpt(k,i) 01430 END DO 01431 END DO 01432 ! 01433 ! Scale D and GD, with a sign change if required. Set S to another 01434 ! vector in the initial two dimensional subspace. 01435 ! 01436 gg=zero 01437 sp=zero 01438 dhd=zero 01439 DO i=1,n 01440 gg=gg+gc(i)**2 01441 sp=sp+d(i)*gc(i) 01442 dhd=dhd+d(i)*gd(i) 01443 END DO 01444 scale=delta/SQRT(dd) 01445 IF (sp*dhd < zero) scale=-scale 01446 temp=zero 01447 IF (sp*sp > 0.99_dp*dd*gg) temp=one 01448 tau=scale*(ABS(sp)+half*scale*ABS(dhd)) 01449 IF (gg*delsq < 0.01_dp*tau*tau) temp=one 01450 DO i=1,n 01451 d(i)=scale*d(i) 01452 gd(i)=scale*gd(i) 01453 s(i)=gc(i)+temp*gd(i) 01454 END DO 01455 ! 01456 ! Begin the iteration by overwriting S with a vector that has the 01457 ! required length and direction, except that termination occurs if 01458 ! the given D and S are nearly parallel. 01459 ! 01460 mainloop : DO 01461 iterc=iterc+1 01462 dd=zero 01463 sp=zero 01464 ss=zero 01465 DO i=1,n 01466 dd=dd+d(i)**2 01467 sp=sp+d(i)*s(i) 01468 ss=ss+s(i)**2 01469 END DO 01470 temp=dd*ss-sp*sp 01471 IF (temp <= 1.0e-8_dp*dd*ss) EXIT mainloop 01472 denom=SQRT(temp) 01473 DO i=1,n 01474 s(i)=(dd*s(i)-sp*d(i))/denom 01475 w(i)=zero 01476 END DO 01477 ! 01478 ! Calculate the coefficients of the objective function on the circle 01479 ! beginning with the multiplication of S by the second derivative ma 01480 ! 01481 DO k=1,npt 01482 sum=zero 01483 DO j=1,n 01484 sum=sum+xpt(k,j)*s(j) 01485 END DO 01486 sum=hcol(k)*sum 01487 DO i=1,n 01488 w(i)=w(i)+sum*xpt(k,i) 01489 END DO 01490 END DO 01491 cf1=zero 01492 cf2=zero 01493 cf3=zero 01494 cf4=zero 01495 cf5=zero 01496 DO i=1,n 01497 cf1=cf1+s(i)*w(i) 01498 cf2=cf2+d(i)*gc(i) 01499 cf3=cf3+s(i)*gc(i) 01500 cf4=cf4+d(i)*gd(i) 01501 cf5=cf5+s(i)*gd(i) 01502 END DO 01503 cf1=half*cf1 01504 cf4=half*cf4-cf1 01505 ! 01506 ! Seek the value of the angle that maximizes the modulus of TAU. 01507 ! 01508 taubeg=cf1+cf2+cf4 01509 taumax=taubeg 01510 tauold=taubeg 01511 isave=0 01512 iu=49 01513 temp=twopi/REAL(iu+1,DP) 01514 DO i=1,iu 01515 angle=REAL(i,dp)*temp 01516 cth=COS(angle) 01517 sth=SIN(angle) 01518 tau=cf1+(cf2+cf4*cth)*cth+(cf3+cf5*cth)*sth 01519 IF (ABS(tau) > ABS(taumax)) THEN 01520 taumax=tau 01521 isave=i 01522 tempa=tauold 01523 ELSE IF (i == isave+1) THEN 01524 tempb=taU 01525 END IF 01526 tauold=tau 01527 END DO 01528 IF (isave == 0) tempa=tau 01529 IF (isave == iu) tempb=taubeg 01530 step=zero 01531 IF (tempa /= tempb) THEN 01532 tempa=tempa-taumax 01533 tempb=tempb-taumax 01534 step=half*(tempa-tempb)/(tempa+tempb) 01535 END IF 01536 angle=temp*(REAL(isave,DP)+step) 01537 ! 01538 ! Calculate the new D and GD. Then test for convergence. 01539 ! 01540 cth=COS(angle) 01541 sth=SIN(angle) 01542 tau=cf1+(cf2+cf4*cth)*cth+(cf3+cf5*cth)*sth 01543 DO i=1,n 01544 d(i)=cth*d(i)+sth*s(i) 01545 gd(i)=cth*gd(i)+sth*w(i) 01546 s(i)=gc(i)+gd(i) 01547 END DO 01548 IF (ABS(tau) <= 1.1_dp*ABS(taubeg)) EXIT mainloop 01549 IF (iterc >= n) EXIT mainloop 01550 END DO mainloop 01551 01552 END SUBROUTINE biglag 01553 !****************************************************************************** 01554 01555 SUBROUTINE trsapp (n,npt,xopt,xpt,gq,hq,pq,delta,step,d,g,hd,hs,crvmin) 01556 01557 INTEGER, INTENT(INOUT) :: n, npt 01558 REAL(dp), DIMENSION(*), INTENT(INOUT) :: xopt 01559 REAL(dp), DIMENSION(npt, *), 01560 INTENT(INOUT) :: xpt 01561 REAL(dp), DIMENSION(*), INTENT(INOUT) :: gq, hq, pq 01562 REAL(dp), INTENT(INOUT) :: delta 01563 REAL(dp), DIMENSION(*), INTENT(INOUT) :: step, d, g, hd, hs 01564 REAL(dp), INTENT(INOUT) :: crvmin 01565 01566 REAL(dp), PARAMETER :: half = 0.5_dp, zero = 0.0_dp 01567 01568 INTEGER :: i, isave, iterc, itermax, 01569 itersw, iu, j 01570 LOGICAL :: jump1, jump2 01571 REAL(dp) :: alpha, angle, angtest, bstep, cf, cth, dd, delsq, dg, dhd, 01572 dhs, ds, gg, ggbeg, ggsav, qadd, qbeg, qmin, qnew, qred, qsav, ratio, 01573 reduc, sg, sgk, shs, ss, sth, temp, tempa, tempb 01574 01575 ! 01576 ! N is the number of variables of a quadratic objective function, Q 01577 ! The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meani 01578 ! in order to define the current quadratic model Q. 01579 ! DELTA is the trust region radius, and has to be positive. 01580 ! STEP will be set to the calculated trial step. 01581 ! The arrays D, G, HD and HS will be used for working space. 01582 ! CRVMIN will be set to the least curvature of H along the conjugate 01583 ! directions that occur, except that it is set to zero if STEP goe 01584 ! all the way to the trust region boundary. 01585 ! 01586 ! The calculation of STEP begins with the truncated conjugate gradient 01587 ! method. If the boundary of the trust region is reached, then further 01588 ! changes to STEP may be made, each one being in the 2D space spanned 01589 ! by the current STEP and the corresponding gradient of Q. Thus STEP 01590 ! should provide a substantial reduction to Q within the trust region 01591 ! 01592 ! Initialization, which includes setting HD to H times XOPT. 01593 ! 01594 01595 delsq=delta*delta 01596 iterc=0 01597 itermax=n 01598 itersw=itermax 01599 DO i=1,n 01600 d(i)=xopt(i) 01601 END DO 01602 CALL updatehd 01603 ! 01604 ! Prepare for the first line search. 01605 ! 01606 qred=zero 01607 dd=zero 01608 DO i=1,n 01609 step(i)=zero 01610 hs(i)=zero 01611 g(i)=gq(i)+hd(i) 01612 d(i)=-g(i) 01613 dd=dd+d(i)**2 01614 END DO 01615 crvmin=zero 01616 IF (dd == zero) RETURN 01617 ds=zero 01618 ss=zero 01619 gg=dd 01620 ggbeg=gg 01621 ! 01622 ! Calculate the step to the trust region boundary and the product HD 01623 ! 01624 jump1 = .FALSE. 01625 jump2 = .FALSE. 01626 mainloop : DO 01627 IF ( .NOT. jump2 ) THEN 01628 IF ( .NOT. jump1 ) THEN 01629 iterc=iterc+1 01630 temp=delsq-ss 01631 bstep=temp/(ds+SQRT(ds*ds+dd*temp)) 01632 CALL updatehd 01633 END IF 01634 jump1 = .FALSE. 01635 IF (iterc <= itersw) THEN 01636 dhd=zero 01637 DO j=1,n 01638 dhd=dhd+d(j)*hd(j) 01639 END DO 01640 ! 01641 ! Update CRVMIN and set the step-length ALPHA. 01642 ! 01643 alpha=bstep 01644 IF (dhd > zero) THEN 01645 temp=dhd/dd 01646 IF (iterc == 1) crvmin=temp 01647 crvmin=MIN(crvmin,temp) 01648 alpha=MIN(alpha,gg/dhd) 01649 END IF 01650 qadd=alpha*(gg-half*alpha*dhd) 01651 qred=qred+qadd 01652 ! 01653 ! Update STEP and HS. 01654 ! 01655 ggsav=gg 01656 gg=zero 01657 DO i=1,n 01658 step(i)=step(i)+alpha*d(i) 01659 hs(i)=hs(i)+alpha*hd(i) 01660 gg=gg+(g(i)+hs(i))**2 01661 END DO 01662 ! 01663 ! Begin another conjugate direction iteration if required. 01664 ! 01665 IF (alpha < bstep) THEN 01666 IF (qadd <= 0.01_dp*qred) EXIT mainloop 01667 IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop 01668 IF (iterc == itermax) EXIT mainloop 01669 temp=gg/ggsav 01670 dd=zero 01671 ds=zero 01672 ss=zero 01673 DO i=1,n 01674 d(i)=temp*d(i)-g(i)-hs(i) 01675 dd=dd+d(i)**2 01676 ds=ds+d(i)*step(I) 01677 ss=ss+step(i)**2 01678 END DO 01679 IF (ds <= zero) EXIT mainloop 01680 IF (ss < delsq) CYCLE mainloop 01681 END IF 01682 crvmin=zero 01683 itersw=iterc 01684 jump2 = .TRUE. 01685 IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop 01686 ELSE 01687 jump2 = .FALSE. 01688 END IF 01689 END IF 01690 ! 01691 ! Test whether an alternative iteration is required. 01692 ! 01693 !!!! IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop 01694 IF (jump2) THEN 01695 sg=zero 01696 shs=zero 01697 DO i=1,n 01698 sg=sg+step(i)*g(i) 01699 shs=shs+step(i)*hs(i) 01700 END DO 01701 sgk=sg+shs 01702 angtest=sgk/SQRT(gg*delsq) 01703 IF (angtest <= -0.99_dp) EXIT mainloop 01704 ! 01705 ! Begin the alternative iteration by calculating D and HD and some 01706 ! scalar products. 01707 ! 01708 iterc=iterc+1 01709 temp=SQRT(delsq*gg-sgk*sgk) 01710 tempa=delsq/temp 01711 tempb=sgk/temp 01712 DO i=1,n 01713 d(i)=tempa*(g(i)+hs(i))-tempb*step(i) 01714 END DO 01715 CALL updatehd 01716 IF (iterc <= itersw) THEN 01717 jump1 = .TRUE. 01718 CYCLE mainloop 01719 END IF 01720 END IF 01721 dg=zero 01722 dhd=zero 01723 dhs=zero 01724 DO i=1,n 01725 dg=dg+d(i)*g(i) 01726 dhd=dhd+hd(i)*d(i) 01727 dhs=dhs+hd(i)*step(i) 01728 END DO 01729 ! 01730 ! Seek the value of the angle that minimizes Q. 01731 ! 01732 cf=half*(shs-dhd) 01733 qbeg=sg+cf 01734 qsav=qbeg 01735 qmin=qbeg 01736 isave=0 01737 iu=49 01738 temp=twopi/REAL(iu+1,dp) 01739 DO i=1,iu 01740 angle=REAL(i,dp)*temp 01741 cth=COS(angle) 01742 sth=SIN(angle) 01743 qnew=(sg+cf*cth)*cth+(dg+dhs*cth)*sth 01744 IF (qnew < qmin) THEN 01745 qmin=qnew 01746 isave=i 01747 tempa=qsav 01748 ELSE IF (i == isave+1) THEN 01749 tempb=qnew 01750 END IF 01751 qsav=qnew 01752 END DO 01753 IF (isave == zero) tempa=qnew 01754 IF (isave == iu) tempb=qbeg 01755 angle=zero 01756 IF (tempa /= tempb) THEN 01757 tempa=tempa-qmin 01758 tempb=tempb-qmin 01759 angle=half*(tempa-tempb)/(tempa+tempb) 01760 END IF 01761 angle=temp*(REAL(isave,DP)+angle) 01762 ! 01763 ! Calculate the new STEP and HS. Then test for convergence. 01764 ! 01765 cth=COS(angle) 01766 sth=SIN(angle) 01767 reduc=qbeg-(sg+cf*cth)*cth-(dg+dhs*cth)*sth 01768 gg=zero 01769 DO i=1,n 01770 step(i)=cth*step(i)+sth*d(i) 01771 hs(i)=cth*hs(i)+sth*hd(i) 01772 gg=gg+(g(i)+hs(i))**2 01773 END DO 01774 qred=qred+reduc 01775 ratio=reduc/qred 01776 IF (iterc < itermax .AND. ratio > 0.01_dp) THEN 01777 jump2 = .TRUE. 01778 ELSE 01779 EXIT mainloop 01780 END IF 01781 01782 IF (gg <= 1.0e-4_dp*ggbeg) EXIT mainloop 01783 01784 END DO mainloop 01785 01786 !******************************************************************************* 01787 CONTAINS 01788 SUBROUTINE updatehd 01789 INTEGER :: i, ih, j, k 01790 01791 DO i=1,n 01792 hd(i)=zero 01793 END DO 01794 DO k=1,npt 01795 temp=zero 01796 DO j=1,n 01797 temp=temp+xpt(k,j)*d(j) 01798 END DO 01799 temp=temp*pq(k) 01800 DO i=1,n 01801 hd(i)=hd(i)+temp*xpt(k,i) 01802 END DO 01803 END DO 01804 ih=0 01805 DO j=1,n 01806 DO i=1,j 01807 ih=ih+1 01808 IF (i < j) hd(j)=hd(j)+hq(ih)*d(i) 01809 hd(i)=hd(i)+hq(ih)*d(j) 01810 END DO 01811 END DO 01812 END SUBROUTINE updatehd 01813 01814 END SUBROUTINE trsapp 01815 !****************************************************************************** 01816 01817 SUBROUTINE update (n,npt,bmat,zmat,idz,ndim,vlag,beta,knew,w) 01818 01819 INTEGER, INTENT(IN) :: n, npt, ndim 01820 INTEGER, INTENT(INOUT) :: idz 01821 REAL(dp), DIMENSION(npt, *), 01822 INTENT(INOUT) :: zmat 01823 REAL(dp), DIMENSION(ndim, *), 01824 INTENT(INOUT) :: bmat 01825 REAL(dp), DIMENSION(*), INTENT(INOUT) :: vlag 01826 REAL(dp), INTENT(INOUT) :: beta 01827 INTEGER, INTENT(INOUT) :: knew 01828 REAL(dp), DIMENSION(*), INTENT(INOUT) :: w 01829 01830 REAL(dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp 01831 01832 INTEGER :: i, iflag, j, ja, jb, jl, jp, 01833 nptm 01834 REAL(dp) :: alpha, denom, scala, scalb, 01835 tau, tausq, temp, tempa, tempb 01836 01837 ! The arrays BMAT and ZMAT with IDZ are updated, in order to shift the 01838 ! interpolation point that has index KNEW. On entry, VLAG contains the 01839 ! components of the vector Theta*Wcheck+e_b of the updating formula 01840 ! (6.11), and BETA holds the value of the parameter that has this na 01841 ! The vector W is used for working space. 01842 ! 01843 01844 nptm=npt-n-1 01845 ! 01846 ! Apply the rotations that put zeros in the KNEW-th row of ZMAT. 01847 ! 01848 jl=1 01849 DO j=2,nptm 01850 IF (j == idz) THEN 01851 jl=idz 01852 ELSE IF (zmat(knew,j) /= zero) THEN 01853 temp=SQRT(zmat(knew,jl)**2+zmat(knew,j)**2) 01854 tempa=zmat(knew,jl)/temp 01855 tempb=zmat(knew,j)/temp 01856 DO I=1,NPT 01857 temp=tempa*zmat(i,jl)+tempb*zmat(i,j) 01858 zmat(i,j)=tempa*zmat(i,j)-tempb*zmat(i,jl) 01859 zmat(i,jl)=temp 01860 END DO 01861 zmat(knew,j)=zero 01862 END IF 01863 END DO 01864 ! 01865 ! Put the first NPT components of the KNEW-th column of HLAG into W, 01866 ! and calculate the parameters of the updating formula. 01867 ! 01868 tempa=zmat(knew,1) 01869 IF (idz >= 2) tempa=-tempa 01870 IF (jl > 1) tempb=zmat(knew,jl) 01871 DO i=1,npt 01872 w(i)=tempa*zmat(i,1) 01873 IF (jl > 1) w(i)=w(i)+tempb*zmat(i,jl) 01874 END DO 01875 alpha=w(knew) 01876 tau=vlag(knew) 01877 tausq=tau*tau 01878 denom=alpha*beta+tausq 01879 vlag(knew)=vlag(knew)-one 01880 ! 01881 ! Complete the updating of ZMAT when there is only one nonzero eleme 01882 ! in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to 01883 ! then the first column of ZMAT will be exchanged with another one l 01884 ! 01885 iflag=0 01886 IF (JL == 1) THEN 01887 temp=SQRT(ABS(denom)) 01888 tempb=tempa/temp 01889 tempa=tau/temp 01890 DO i=1,npt 01891 zmat(i,1)=tempa*zmat(i,1)-tempb*vlag(i) 01892 END DO 01893 IF (idz == 1 .AND. temp < zero) idz=2 01894 IF (idz >= 2 .AND. temp >= zero) iflag=1 01895 ELSE 01896 ! 01897 ! Complete the updating of ZMAT in the alternative case. 01898 ! 01899 ja=1 01900 IF (beta >= zero) ja=jl 01901 jb=jl+1-ja 01902 temp=zmat(knew,jb)/denom 01903 tempa=temp*beta 01904 tempb=temp*tau 01905 temp=zmat(knew,ja) 01906 scala=one/SQRT(ABS(beta)*temp*temp+tausq) 01907 scalb=scala*SQRT(ABS(denom)) 01908 DO i=1,npt 01909 zmat(i,ja)=scala*(tau*zmat(i,ja)-temp*vlag(i)) 01910 zmat(i,jb)=scalb*(zmat(i,jb)-tempa*w(i)-tempb*vlag(i)) 01911 END DO 01912 IF (denom <= zero) THEN 01913 IF (beta < zero) idz=idz+1 01914 IF (beta >= zero) iflag=1 01915 END IF 01916 END IF 01917 ! 01918 ! IDZ is reduced in the following case, and usually the first column 01919 ! of ZMAT is exchanged with a later one. 01920 ! 01921 IF (iflag == 1) THEN 01922 idz=idz-1 01923 DO i=1,npt 01924 temp=zmat(i,1) 01925 zmat(i,1)=zmat(i,idz) 01926 zmat(i,idz)=temp 01927 END DO 01928 END IF 01929 ! 01930 ! Finally, update the matrix BMAT. 01931 ! 01932 DO j=1,n 01933 jp=npt+j 01934 w(jp)=bmat(knew,j) 01935 tempa=(alpha*vlag(jp)-tau*w(jp))/denom 01936 tempb=(-beta*w(jp)-tau*vlag(jp))/denom 01937 DO i=1,jp 01938 bmat(i,j)=bmat(i,j)+tempa*vlag(i)+tempb*w(i) 01939 IF (i > npt) bmat(jp,i-npt)=bmat(i,j) 01940 END DO 01941 END DO 01942 01943 END SUBROUTINE update 01944 01945 END MODULE powell 01946 01947 !******************************************************************************
1.7.3