C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       LMINLS        **
C*************************
C
C
      SUBROUTINE LMINLS(P,XI,N,FRET,M,FVECC,XTEMP,CALFUN,IFAIL)
C
C Purpose: Line mimimisation routine for use with least squares mimimisation
C          routine SSQMIN
C
C          Based on methods and code described in 'Numerical Recipes in C'
C          by W. Press, B. Flannery, S. Teukolsky and W. Vettering
C          Cambridge University Press
C
      DOUBLE PRECISION P(N), XI(N), FVECC(M), XTEMP(N)
      DOUBLE PRECISION FRET
      INTEGER N, M
      EXTERNAL CALFUN
      INTEGER IFAIL

      DOUBLE PRECISION FNCALC
C
C====== LOCALS:
C
      INTEGER MAXITR

      DOUBLE PRECISION AX, BX, CX
      DOUBLE PRECISION FA, FB, FC, TEMP
      DOUBLE PRECISION Q, R, U, ULIM, FU
      DOUBLE PRECISION A, B, E, D, ETEMP
      DOUBLE PRECISION X, XMIN, XMEAN
      DOUBLE PRECISION PP, QQ, RR, V, W
      DOUBLE PRECISION FX, FV, FW
      DOUBLE PRECISION TOLER1, TOLER2

      DOUBLE PRECISION GOLD, CGOLD, GLIMIT, TINY, ZEPS, TOLER
      DATA GOLD/1.618034/
      DATA CGOLD/.381966/
      DATA GLIMIT/100.0/
      DATA TINY/1.0D-20/
      DATA ZEPS/1.0D-10/
      DATA TOLER/2.0D-4/
C
C-------------------------------------------------------------------------------
C
C
C====== Initialisations
C
      MAXITR = 100
      IFAIL = 0
      AX=0.0
      BX=1.0
      CX=2.0


C
C====== Get bracketing values
C


C
C====== Calculate the function offset at AX and BX
C
      FA=FNCALC(P,XI,AX,M,N,FVECC,XTEMP,CALFUN)
      FB=FNCALC(P,XI,BX,M,N,FVECC,XTEMP,CALFUN)
C
C====== Swap points if value at BX greater than value at AX
C
      IF (FB.GT.FA) THEN
         TEMP=AX
         AX=BX
         BX=TEMP
         TEMP=FB
         FB=FA
         FA=TEMP
      END IF
C
C====== Set initial value for CX
C
      CX = BX + GOLD*(BX-AX)
      FC = FNCALC(P,XI,CX,M,N,FVECC,XTEMP,CALFUN)
C
C====== Iteration loop
C
10    IF (FB.LE.FC) GO TO 100
      R=(BX-AX)*(FB-FC)
      Q=(BX-CX)*(FB-FA)
      U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.0*SIGN(MAX(ABS(Q-R),TINY),Q-R))
      ULIM=BX+GLIMIT*(CX-BX)
      IF ((BX-U)*(U-CX).GT.0.0) THEN
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
         IF (FU.LT.FC) THEN
            AX=BX
            BX=U
            FA=FB
            FB=FU
            GO TO 100
         ELSE IF (FU.GT.FB) THEN
            CX=U
            FC=U
            GO TO 100
         END IF
         U=CX+GOLD*(CX-BX)
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
      ELSE IF ((CX-U)*(U-ULIM).GT.0.0) THEN
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
         IF (FU.LT.FC) THEN
            BX=CX
            CX=U
            U=CX+GOLD*(CX-BX)
            FB=FC
            FC=FU
            FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
         END IF
      ELSE IF ((U-ULIM)*(ULIM-CX).GE.0.0) THEN
         U=ULIM
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
      ELSE
         U=CX+GOLD*(CX-BX)
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
      END IF
      AX=BX
      BX=CX
      CX=U
      FA=FB
      FB=FC
      FC=FU
      GO TO 10


C
C====== Find the mimimum using Brent's method
C


100   E=0.0
      D=0.0
C
C====== Get limits a,b in ascending order
C
      IF (AX.LT.CX) THEN
         A=AX
         B=CX
      ELSE
         A=CX
         B=AX
      ENDIF
      X=BX
      W=X
      V=X
C
C====== Calculate function at the intermediate point
C
      FX=FNCALC(P,XI,X,M,N,FVECC,XTEMP,CALFUN)
      FW=FX
      FV=FX   
C
C====== Iteration Loop
C
      DO 200 ITER=1,MAXITR
         
         XMEAN=0.5*(A+B)
         TOLER1=TOLER*ABS(X)+ZEPS
         TOLER2=2.0*TOLER1
C
C====== See if answer reached
C
         IF (ABS(X-XMEAN).LE.TOLER2-0.5*(B-A)) THEN
            XMIN=X
            FRET=FX
            GO TO 250
         END IF
         IF (ABS(E).GT.TOLER1) THEN
C
C====== Try a parabolic fit
C
            RR=(X-W)*(FX-FV)
            QQ=(X-V)*(FX-FW)
            PP=(X-V)*QQ-(X-W)*RR
            QQ=2.0*(QQ-RR)
            IF (QQ.GT.0.0) PP=-PP
            QQ=ABS(QQ)
            ETEMP=E
            E=D
            IF (ABS(PP).GE.ABS(0.5*QQ*ETEMP)
     +          .OR.PP.LE.QQ*(A-X)
     +          .OR.PP.GE.QQ*(B-X)) THEN
C
C====== Take the golden section step into the larger of the two segments
C
               IF (X.GE.XMEAN) THEN
                  E=A-X
               ELSE
                  E=B-X
               END IF
               D=CGOLD*E
            ELSE
C
C====== Take the parabolic step
C
               D=PP/QQ
               U=X+D
               IF (U-A.LT.TOLER2.OR.B-U.LT.TOLER2) THEN
                  D=SIGN(TOLER1,XMEAN-X)
               END IF
            END IF
         ELSE
            IF (X.GE.XMEAN) THEN
               E=A-X
            ELSE
               E=B-X
            ENDIF
            D=CGOLD*E
         END IF
         IF (ABS(D).GE.TOLER1) THEN
            U=X+D
         ELSE
            U=X+SIGN(TOLER1,D)
         END IF
C
C====== Evaluate the function and prepare for next iteration
C
         FU=FNCALC(P,XI,U,M,N,FVECC,XTEMP,CALFUN)
         IF (FU.LE.FX) THEN
            IF (U.GE.X) THEN
               A=X
            ELSE
               B=X
            ENDIF
            V=W
            W=X
            X=U
            FV=FW
            FW=FX
            FX=FU
         ELSE
            IF (U.LT.X) THEN
               A=U
            ELSE
               B=U
            END IF
            IF (FU.LE.FW.OR.W.EQ.X) THEN
               V=W
               W=U
               FV=FW
               FW=FU
            ELSE IF (FU.LE.FV.OR.V.EQ.X.OR.V.EQ.W) THEN
               V=U
               FV=FU
            END IF
         END IF
200   CONTINUE
      IFAIL=1
      XMIN=X
      FRET=FX
C
C====== Return vector of results
C
250   DO 260 I=1,N
         XI(I)=XI(I)*XMIN
         P(I)=P(I)+XI(I)
260   CONTINUE
      RETURN
      END