C
      SUBROUTINE POWLL(FUNC,XVLO,XVEC,XVHI,NVEC)
C
C Subroutine POWLL: Powell-method minimization: uses the
C                   NR subroutines of Press et al, with a
C                   slight style modification.
C
C Update history:
C    V1.0    Original version                            8-Oct-1992
C
      IMPLICIT NONE
C
C Arguments:
C    FUNC = R*4 function = quantity to minimize
C    XVLO = R*4(1)       = array of lower limits to XVEC
C    XVEC = R*4(1)       = array of quantities to find at minimum
C    XVHI = R*4(1)       = array of upper limits to XVEC
C    NVEC = I*4          = number of parameters
C
      REAL*4    FUNC
      EXTERNAL  FUNC
C
      REAL*4    XVLO(1),XVEC(1),XVHI(1)
      INTEGER*4 NVEC
C
C Subprogrammes:
C    LINMIN = subroutine = modified Press et al routine
C
      EXTERNAL LINMIN
C
C Local parameter:
C    NVMAX = I*4 = maximum length of vectors
C
      INTEGER*4  NVMAX
      PARAMETER (NVMAX=    128)
C
C Local variables:
C    DEL       = R*4              = amount of function change
C    EPSPW     = R*4              = fractional accuracy demanded
C    FMIN      = R*4              = minimum of function
C    FP        = R*4              = value of function at current point
C    FPTT      = R*4              = value at previous point
C    I1        = I*4              = counting index
C    I2        = I*4              = counting index
C    IBIG      = I*4              = direction of largest change
C    ITER      = I*4              = number of iterations used
C    NIMAX     = I*4              = maximum allowed iterations
C    TEST      = R*4              = test quantity in method
C    PTT(i)    = R*4(NVMAX)       = previous parameters
C    XINI(i)   = R*4(NVMAX)       = initial search direction
C    XIT(i)    = R*4(NVMAX)       = amount of shift
C    XMAT(i,j) = R*4(NVMAX,NVMAX) = matrix used by POWELL
C
      REAL*4    DEL,EPSPW,FMIN,FP,FPTT,TEST
      REAL*4    PTT(NVMAX),XINI(NVMAX),XIT(NVMAX),XMAT(NVMAX,NVMAX)
      INTEGER*4 I1,I2,IBIG,ITER,NIMAX
C
C Check that vector length isn't too large
C
      IF (NVEC .GT. NVMAX) THEN
        WRITE (6,699)
        RETURN
      ENDIF
C
C Get controlling variables
C
      CALL GETOP('POWLL.EPS',EPSPW)
      CALL GETOP('POWLL.NIMAX',TEST)
      NIMAX=IFIX(TEST+0.5)
C
C Initialise matrix of search directions
C
      DO I1=1,NVEC
        DO I2=1,NVEC
          XMAT(I2,I1)=0.0
        ENDDO
        XMAT(I1,I1)=1.0
      ENDDO
C
C Use initial vector to get initial minimum; save initial vector
C
      FMIN=FUNC(XVEC,NVEC)
      DO I1=1,NVEC
        XINI(I1)=XVEC(I1)
      ENDDO
C
      WRITE (6,999) FMIN
C
C Initialise counter of interations
C
      ITER=0
C
C Main loop, done once per iteration
C
 1    CONTINUE
        ITER=ITER+1
C
C Save current value of the minimum, and initialise the pointer
C    to the `best' direction in which to do the minimum (IBIG).
C    Also initialise DEL, the largest value of the decrease in
C    the function.
C
        FP  =FMIN
        IBIG=0
        DEL =0.0
C
C Loop over directions in the matrix, stepping once in each of
C    the directions to locate the minimum.
C
        DO I1=1,NVEC
C
C Select direction into vector XIT
C
          DO I2=1,NVEC
            XIT(I2)=XMAT(I2,I1)
          ENDDO
C
C Find the minimum in the function in direction XIT, allowing 
C    motion only within the legal bounds.
C
          CALL LINMIN(FUNC,XVLO,XVEC,XVHI,XIT,NVEC,FMIN)
C
C Record the direction index which gives the largest
C    decrease in the function
C
          IF (ABS(FP-FMIN) .GT. DEL) THEN
            DEL =ABS(FP-FMIN)
            IBIG=I1
          ENDIF
        ENDDO
C
C Now, if the function decrease is less than the
C    tolerance, the job is done.
C
        IF (ABS(FP-FMIN) .LE. 0.5*EPSPW*(ABS(FP)+ABS(FMIN))) THEN
          WRITE (6,998) FMIN,ITER
          RETURN
        ENDIF
C
C If the number of iterations has reached the limit, then
C    the job is done
C
        IF (ITER .GE. NIMAX) THEN
          WRITE (6,997) FMIN,ITER
          RETURN
        ENDIF
C
C Construct the extrapolated point to which the minimum moved,
C    and store as the initial seach point for the next loop.
C    
        DO I2=1,NVEC
          PTT(I2) =AMAX0(XVLO(I2),
     +                   AMIN0(2.0*XVEC(I2)-XINI(I2),XVHI(I2)))
          XIT(I2) =    XVEC(I2)-XINI(I2)
          XINI(I2)=    XVEC(I2)
        ENDDO
C
C Get value of function at extrapolated point beyond minimum
C
        FPTT=FUNC(PTT,NVEC)
C
C If this value is greater than the original value at the start
C    of the loop, then don't include the new direction in the
C    search
C
        IF (FPTT .GE. FP) THEN
          WRITE (6,996) FMIN,ITER
          GOTO 1
        ENDIF
C
C If near bottom of minimum, or the average direction is not 
C    particularly good, don't include new direction in search
C
	TEST=2.0*(FP-2.0*FMIN+FPTT)*(FP-FMIN-DEL)**2
     +       -DEL*(FP-FPTT)**2
	IF (TEST .GE. 0.0) THEN
          WRITE (6,995) FMIN,ITER
          GOTO 1
        ENDIF
C
C Locate the minimum in the new direction created from the
C    search
C
        CALL LINMIN(FUNC,XVLO,XVEC,XVHI,XIT,NVEC,FMIN)
C
C Save the new direction into the search matrix to replace the
C    previous best direction
C
        DO I2=1,NVEC
          XMAT(I2,IBIG)=XIT(I2)
	ENDDO
C
      WRITE (6,994) FMIN,ITER
C
C Return to top of loop
C
      GOTO 1
C
C Error messages
C
 699  FORMAT (' '/
     +        ' *** Error: parameter vector too large ***'/)
C
C Debug progress messages
C
 994  FORMAT (' POWLL: (new direction) at value = ',
     +        1PE15.5,' at iteration = ',I4)
 995  FORMAT (' POWLL:   (failed TEST) at value = ',
     +        1PE15.5,' at iteration = ',I4)
 996  FORMAT (' POWLL:   (failed FPTT) at value = ',
     +        1PE15.5,' at iteration = ',I4)
 997  FORMAT (' POWLL: iteration limit at value = ',
     +        1PE15.5,' at iteration = ',I4)
 998  FORMAT (' POWLL:     converged to minimum = ',
     +        1PE15.5,' at iteration = ',I4)
 999  FORMAT (' '/
     +        ' POWLL:   initial function value = ',
     +        1PE15.5)
C
      END