C ***************************************************************** C ***************************************************************** subroutine extrapolation(nind,x,l,u,m,lambda,rho,equatn,linear,g, +xp,fp,gp,d,alpha,amax,rbdnnz,rbdind,rbdtype,fmin,beta,etaext, +maxextrap,extinfo,inform) C SCALAR ARGUMENTS integer extinfo,inform,m,maxextrap,nind,rbdnnz double precision alpha,amax,beta,etaext,fmin,fp C ARRAY ARGUMENTS logical equatn(m),linear(m) character rbdtype(nind) integer rbdind(nind) double precision d(nind),g(*),gp(*),l(nind),lambda(m),rho(m), + u(nind),x(*),xp(*) C Performs extrapolation. C C extinfo: C C 0: Success C 2: Unbounded objective function? C 4: beta-condition holds. No extrapolation is done C 6: Maximum number of extrapolations reached C 7: Similar consecutive projected points C 8: Not-well-defined objective function C 9: Functional value increases include "dim.par" include "counters.inc" include "machconst.inc" include "outtyp.inc" C LOCAL SCALARS logical projected,samep integer extrap,i double precision atmp,fbext,ftmp,gptd,gtd C LOCAL ARRAYS double precision xbext(nmax),xtmp(nmax) extinfo = 0 extrap = 0 C Compute directional derivative gtd = 0.0d0 do i = 1,nind gtd = gtd + g(i) * d(i) end do call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return gptd = 0.0d0 do i = 1,nind gptd = gptd + gp(i) * d(i) end do C Beta condition holds. No extrapolation is done. if ( gptd .ge. beta * gtd ) then if ( iprintinn .ge. 6 ) then C write(*, 110) write(10,110) end if extinfo = 4 return end if C Beta condition does not holds. We will extrapolate. if ( iprintinn .ge. 6 ) then C write(* ,120) write(10,120) end if C Save f and x before extrapolation to return in case of a C not-well-defined objective function at an extrapolated point fbext = fp do i = 1,nind xbext(i) = xp(i) end do 1010 continue C Test f going to -inf if ( fp .le. fmin ) then C Finish the extrapolation with the current point if ( extrap .ne. 0 ) then call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return end if extinfo = 2 if ( iprintinn .ge. 6 ) then C write(*, 910) write(10,910) end if return end if C Test if the maximum number of extrapolations was exceeded if ( extrap .ge. maxextrap ) then C Finish the extrapolation with the current point if ( extrap .ne. 0 ) then call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return end if extinfo = 6 if ( iprintinn .ge. 6 ) then C write(*, 930) write(10,930) end if return end if C Chose new step extrap = extrap + 1 if ( alpha .lt. amax .and. etaext * alpha .gt. amax ) then atmp = amax else atmp = etaext * alpha end if C Compute new trial point do i = 1,nind xtmp(i) = x(i) + atmp * d(i) end do if ( atmp .eq. amax ) then do i = 1,rbdnnz if ( rbdtype(i) .eq. 'L' ) then xtmp(rbdind(i)) = l(rbdind(i)) else if ( rbdtype(i) .eq. 'U' ) then xtmp(rbdind(i)) = u(rbdind(i)) end if end do end if C Project if ( atmp .gt. amax ) then projected = .false. do i = 1,nind if ( xtmp(i) .lt. l(i) .or. xtmp(i) .gt. u(i) ) then projected = .true. xtmp(i) = max( l(i), min( xtmp(i), u(i) ) ) end if end do end if C Test if this is not the same point as the previous one. This test C is performed only when xtmp is in fact a projected point. if ( projected ) then samep = .true. do i = 1,nind if ( xtmp(i) .gt. + xp(i) + macheps23 * max( 1.0d0, abs( xp(i) ) ) .or. + xtmp(i) .lt. + xp(i) - macheps23 * max( 1.0d0, abs( xp(i) ) ) ) then samep = .false. end if end do if ( samep ) then C Finish the extrapolation with the current point call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return extinfo = 7 if ( iprintinn .ge. 6 ) then C write(*, 940) write(10,940) end if return end if end if C Evaluate function call csetp(nind,xtmp) call calcal(nind,xtmp,m,lambda,rho,equatn,linear,ftmp,inform) C If the objective function is not well defined in an extrapolated C point, we discard all the extrapolated points and return to a C safe region (to the last point before the extrapolation) if ( inform .lt. 0 ) then fp = fbext do i = 1,nind xp(i) = xbext(i) end do call csetp(nind,xp) call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return extinfo = 8 if ( iprintinn .ge. 6 ) then C write(*, 950) write(10,950) end if return end if C Print information of this iteration if ( iprintinn .ge. 6 ) then C write(*, 100) atmp,ftmp,fcnt write(10,100) atmp,ftmp,fcnt end if C If the functional value decreases then set the current point and C continue the extrapolation if ( ftmp .lt. fp ) then alpha = atmp fp = ftmp do i = 1,nind xp(i) = xtmp(i) end do go to 1010 end if C If the functional value does not decrease then discard the last C trial and finish the extrapolation with the previous point call csetp(nind,xp) call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return extinfo = 9 if ( iprintinn .ge. 6 ) then C write(*, 960) write(10,960) end if C NON-EXECUTABLE STATEMENTS 100 format( 5X,'Alpha = ',1P,D7.1,' F = ',1P,D24.16,' FE = ',I7) 110 format( 5X,'Beta condition also holds. ', + 'No extrapolation is done.') 120 format( 5X,'Beta condition does not hold. We will extrapolate.') 910 format( 5X,'Flag of Extrapolation: Unbounded objective ', + 'function?') 930 format( 5X,'Flag of Extrapolation: Maximum of consecutive ', + 'extrapolations reached.') 940 format( 5X,'Flag of Extrapolation: Very similar consecutive ', + 'projected points.') 950 format( 5X,'Flag of Extrapolation: Not-well-defined objective ', + 'function in an extrapolated point.') 960 format( 5X,'Flag of Extrapolation: Functional value increased ', + 'when extrapolating.') end