C ***************************************************************** C ***************************************************************** subroutine tnls(nind,x,l,u,m,lambda,rho,equatn,linear,f,g,amax,d, +rbdnnz,rbdind,rbdtype,xp,fp,gp,lsinfo,inform) implicit none C SCALAR ARGUMENTS integer inform,lsinfo,m,nind,rbdnnz double precision amax,f,fp C ARRAY ARGUMENTS integer rbdind(nind) character rbdtype(nind) logical equatn(m),linear(m) double precision d(nind),g(*),gp(*),l(nind),lambda(m),rho(m), + u(nind),x(*),xp(*) C This subroutine computes a line serach in direction d. C C lsinfo: C C At the first trial: C C 4: Interior first trial. Armijo and beta-condition hold. C 5: x + amax d is at the boundary and f(x + amax d) is smaller than f C Extrapolation: C C 2: Unbounded objective function? 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 C Backtracking: C C 0: Armijo satisfied C 1: Small step with functional value similar to the current one C 2: Unbounded objective function? C 3: Too small backtracking step. Wrong gradient? #include "dim.par" #include "machconst.com" #include "algconst.par" #include "counters.com" #include "outtyp.com" C LOCAL SCALARS logical boundary,projected,samep integer extrap,i double precision alpha,atmp,dsupn,fbext,ftmp,gptd,gtd,xsupn C LOCAL ARRAYS double precision xbext(nmax),xtmp(nmax) C EXTERNAL SUBROUTINES external calcal C ================================================================== C ================================================================== C Test Armijo condition and beta condition. Decide between accept C the first trial, extrapolate or to do backtracking. C ================================================================== C ================================================================== C ------------------------------------------------------------------ C Compute directional derivative, dsupn and xsupn C ------------------------------------------------------------------ gtd = 0.0d0 dsupn = 0.0d0 xsupn = 0.0d0 do i = 1,nind gtd = gtd + g(i) * d(i) dsupn = max( dsupn, abs( d(i) ) ) xsupn = max( xsupn, abs( x(i) ) ) end do if ( iprintinn .ge. 6 ) then write(*, 100) xsupn,amax,dsupn write(10,100) xsupn,amax,dsupn end if C ------------------------------------------------------------------ C Compute first trial (projected point) C ------------------------------------------------------------------ alpha = 1.0d0 boundary = .false. do i = 1,nind xp(i) = x(i) + d(i) if ( xp(i) .lt. l(i) .or. xp(i) .gt. u(i) ) then boundary = .true. xp(i) = max( l(i), min( xp(i), u(i) ) ) end if end do if ( amax .eq. 1.0d0 ) then boundary = .true. do i = 1,rbdnnz if ( rbdtype(i) .eq. 'L' ) then xp(rbdind(i)) = l(rbdind(i)) else if ( rbdtype(i) .eq. 'U' ) then xp(rbdind(i)) = u(rbdind(i)) end if end do end if call calcal(nind,xp,m,lambda,rho,equatn,linear,fp,inform) if ( inform .lt. 0 ) return if ( .not. boundary ) then if ( iprintinn .ge. 6 ) then write(*, 110) fp,fcnt write(10,110) fp,fcnt end if else if ( iprintinn .ge. 6 ) then write(*, 120) fp,fcnt write(10,120) fp,fcnt end if end if C ------------------------------------------------------------------ C The first trial is an interior point. C ------------------------------------------------------------------ if ( .not. boundary ) then if ( iprintinn .ge. 6 ) then write(*, 140) write(10,140) end if C Armijo condition holds. if ( fp .le. f + alpha * gamma * gtd ) then if ( iprintinn .ge. 6 ) then write(*, 150) write(10,150) end if 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. Line search is over. if ( gptd .ge. beta * gtd ) then if ( iprintinn .ge. 6 ) then write(*, 160) write(10,160) end if lsinfo = 4 if ( iprintinn .ge. 6 ) then write(*, 900) write(10,900) end if return end if C Beta condtion does not holds. We will extrapolate. if ( iprintinn .ge. 6 ) then write(* ,170) write(10,170) end if go to 1000 end if C Armijo condition does not hold. We will do backtracking. if ( iprintinn .ge. 6 ) then write(* ,180) write(10,180) end if go to 2000 end if C ------------------------------------------------------------------ C The first trial is at the boundary. C ------------------------------------------------------------------ if ( iprintinn .ge. 6 ) then write(*, 190) write(10,190) end if C Function value is smaller than at the current point. We will C extrapolate. if ( fp .lt. f ) then if ( iprintinn .ge. 6 ) then write(*, 200) write(10,200) end if go to 1000 end if C Discard the projected point and consider x + amax d if ( iprintinn .ge. 6 ) then write(*, 210) write(10,210) end if alpha = amax do i = 1,nind xp(i) = x(i) + alpha * d(i) end do do i = 1,rbdnnz if ( rbdtype(i) .eq. 'L' ) then xp(rbdind(i)) = l(rbdind(i)) else if ( rbdtype(i) .eq. 'U' ) then xp(rbdind(i)) = u(rbdind(i)) end if end do call calcal(nind,xp,m,lambda,rho,equatn,linear,fp,inform) if ( inform .lt. 0 ) return if ( iprintinn .ge. 6 ) then write(*, 130) alpha,fp,fcnt write(10,130) alpha,fp,fcnt end if C Function value is smaller than or equal to (or even just a little C bit greater than) at the current point. Line search is over. if ( fp .lt. f + macheps23 * abs(f) ) then if ( iprintinn .ge. 6 ) then write(*, 220) write(10,220) end if call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return lsinfo = 5 if ( iprintinn .ge. 6 ) then write(*, 900) write(10,900) end if return end if C Function value is greater than at the current point. We will C do backtracking. if ( iprintinn .ge. 6 ) then write(*, 230) write(10,230) end if go to 2000 C ================================================================== C ================================================================== C Extrapolation C ================================================================== C ================================================================== 1000 continue extrap = 0 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 call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return lsinfo = 2 if ( iprintinn .ge. 6 ) then 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 call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return lsinfo = 6 if ( iprintinn .ge. 6 ) then 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 ( abs( xtmp(i) - xp(i) ) .gt. + macheps * 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 lsinfo = 7 if ( iprintinn .ge. 6 ) then write(*, 940) write(10,940) end if return end if end if C Evaluate function 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 lsinfo = 8 if ( iprintinn .ge. 6 ) then write(*, 950) write(10,950) end if return end if C Print information of this iteration if ( iprintinn .ge. 6 ) then write(*, 130) atmp,ftmp,fcnt write(10,130) 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 lsinfo = 9 if ( iprintinn .ge. 6 ) then write(*, 960) write(10,960) end if return C ================================================================== C End of extrapolation C ================================================================== C ================================================================== C ================================================================== C Backtracking C ================================================================== C ================================================================== 2000 continue call backtracking(nind,x,m,lambda,rho,equatn,linear,f,xsupn,d,gtd, +dsupn,alpha,fp,xp,calcal,lsinfo,inform) if ( inform .lt. 0 ) return call calcnal(nind,xp,m,lambda,rho,equatn,linear,gp,inform) if ( inform .lt. 0 ) return C ================================================================== C End of backtracking C ================================================================== C NON-EXECUTABLE STATEMENTS 100 format(/,5X,'TN Line search (xsupn = ',1P,D7.1,', amax = ', + 1P,D7.1,', dsupn = ',1P,D7.1,')') 110 format( 5X,'Unitary step F = ',1P,D24.16,' FE = ',I7) 120 format( 5X,'Projected point F = ',1P,D24.16,' FE = ',I7) 130 format( 5X,'Alpha = ',1P,D7.1,' F = ',1P,D24.16,' FE = ',I7) 140 format( 5X,'The first trial is an interior point.') 150 format( 5X,'Armijo condition holds.') 160 format( 5X,'Beta condition also holds. Line search is over.') 170 format( 5X,'Beta condition does not hold. We will extrapolate.') 180 format( 5X,'Armijo condition does not hold. We will backtrack.') 190 format( 5X,'The first trial is at the boundary.') 200 format( 5X,'Function value at the boundary is smaller than at ', + 'the current point.',/,5X,'We will extrapolate.') 210 format( 5X,'Discarding projected point. We will now consider x ', + '+ amax d.') 220 format( 5X,'Function value at the boundary is smaller than or ', + 'equal to than at the',/,5X,'current point. Line ', + 'search is over.') 230 format( 5X,'Function value at the boundary is greater than at ', + 'the current point.',/,5X,'We will backtrack.') 900 format( 5X,'Flag of TN Line search: First trial accepted.') 910 format( 5X,'Flag of TN Line search: Unbounded objective ', + 'function?') 930 format( 5X,'Flag of TN Line search: Maximum of consecutive ', + 'extrapolations reached.') 940 format( 5X,'Flag of TN Line search: Very similar consecutive ', + 'projected points.') 950 format( 5X,'Flag of TN Line search: Not-well-defined objective ', + 'function in an extrapolated point.') 960 format( 5X,'Flag of TN Line search: Functional value increased ', + 'when extrapolating.') end