C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C*************************
C**       SSQMIN        **
C*************************
C
C
      SUBROUTINE SSQMIN(M,N,X,ACC,MAXITR,CALFUN,FMIN,FVECC,IFLAG,W)
C
C Purpose: Minimise the sum of squares of 'm' given functions of 'n'
C          parameters
C
C          The function minimised is:
C
C                      M
C             F(x) = Sigma [f(j)(X(1),X(2),...X(N))]**2         (M>=N)
C                     j=1
C
C          The user must supply an initial approximation to the required
C          solution in X and a subroutine CALFUN to calculate the values
C          of the functions f(j), j=1,2,...M for any values of the parameters
C          in X. (see below for parameters of CALFUN routine)
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 for a modified Powell minimisation
C
C Parameters:
C
C INTEGER               M (I)    The number of functions f(j) occuring in the 
C                                sum of squares.
C INTEGER               N (I)    The number of variables X(1) ... X(N)
C DOUBLE PRECISION      X (I/O)  An array of at least length N containing,
C                                on input, the user's estimate of the
C                                required parameters X(1) ... X(N)
C                                On output the refined parameter values
C                                are returned.
C DOUBLE PRECISION    ACC (I)    Fractional tolerance in the Function value
C                                such that failure to decrease by more than
C                                this amount one one iteration signals
C                                completness.
C INTEGER          MAXITR (I)    The maximum number of iterations to be taken
C EXTERNAL         CALFUN (I)    Subroutine for calculating the function
C                                values for a given set of parameters (see
C                                below)
C DOUBLE PRECISION   FMIN (O)    The returned minimum value sum of squares
C DOUBLE PRECISION  FVECC (O)    An array, dimensioned to at least M, returning
C                                the M values of f(j), j=1...M calculated using
C                                the refined parameter values.
C INTEGER           IFLAG (O)    If >0 returns the number of iterations done.
C                                   -1 one or more failures to find a line
C                                      minimum.
C                                   -2 refinement not yet completed after
C                                      MAXITR iterations.
C DOUBLE PRECISION      W (*)    An array for work space dimensioned to at
C                                least N*N + 4*N.
C
C     The user written function CALFUN has the call
C
C                      CALFUN (M,N,FVECC,X)
C
C          Where the parameters have the same meaning as for SSQMIN; this
C          routine must not alter the values of M, N, or X.
C  
      INTEGER M, N, MAXITR, IFLAG
      DOUBLE PRECISION FVECC(M), X(N)
      DOUBLE PRECISION ACC, FMIN
      EXTERNAL CALFUN

      DOUBLE PRECISION W(N*N+4*N)
C
C Work array split up as follows:
C
C     XI(N,N)   starts at W(IXI+1)
C     PT(N)     starts at W(IPT+1)
C     PTT(N)    starts at W(IPTT+1)
C     XIT(N)    starts at W(IXIT+1)
C     XTEMP(N)  starts at W(IXTEMP+1)
C
C====== LOCALS:
C
      INTEGER I, J, ITER, IFAIL
      INTEGER IXI, IPT, IPTT, IXIT, IXTEMP
      INTEGER IBIG, JFAIL
      DOUBLE PRECISION FP, FPTT, DEL, T
C
C-------------------------------------------------------------------------------
C
C
C====== Set offsets in work array
C
      IXI=0
      IPT=IXI+N*N
      IPTT=IPT+N
      IXIT=IPTT+N
      IXTEMP=IXIT+N
C
C====== Initialise
C
      IFAIL=0
      DO 10 I=1,N
         W(IPT+I)=X(I)
10    CONTINUE
      DO 20 I=1,N*N
         W(IXI+I)=0.0
20    CONTINUE 
      DO 30 I=1,N
         W(IXI+N*(I-1)+I)=1.0
30    CONTINUE
      CALL CALFUN(M,N,FVECC,X)
      FMIN=0.0
      DO 40 I=1,M
         FMIN=FMIN+FVECC(I)*FVECC(I)
40    CONTINUE
C
C====== Start of iterations loop
C
      DO 500 ITER=1,MAXITR
         FP=FMIN
         IBIG=0
         DEL=0.0
         DO 200 I=1,N
            DO 120 J=1,N
               W(IXIT+J)=W(IXI+N*(J-1)+I)
120         CONTINUE
            FPTT=FMIN
            CALL LMINLS(X,W(IXIT+1),N,FMIN,M,FVECC,W(IXTEMP+1),CALFUN,
     +                  JFAIL)
            IF (JFAIL.EQ.1) IFAIL=1
            IF (ABS(FPTT-FMIN).GT.DEL) THEN
               DEL=ABS(FPTT-FMIN)
               IBIG=I
            ENDIF
200      CONTINUE
C
C====== See if termination criterion met
C
         IF (2.0*ABS(FP-FMIN).LE.ACC*(ABS(FP)+ABS(FMIN))) THEN
            IFLAG=ITER
            IF (IFAIL.EQ.1) IFLAG=-1
            RETURN
         ENDIF
C
C======= Calculate extrapolated point and average direction moved
C======= and save the old starting point
C
         DO 220 J=1,N
            W(IPTT+J)=2.0*X(J)-W(IPT+J)
            W(IXIT+J)=X(J)-W(IPT+J)
            W(IPT+J)=X(J)
220      CONTINUE
C
C====== Calculate function at extrapolated point
C
         CALL CALFUN(M,N,FVECC,W(IPTT+1))
         FPTT = 0.0
         DO 230 J=1,M
            FPTT=FPTT+FVECC(J)*FVECC(J)
230      CONTINUE
C
C====== Move as required
C
         IF (FPTT.LT.FP) THEN
            T=2.0*(FMIN+FPTT)*((FP-FMIN-DEL)**2)-DEL*((FP-FPTT)**2)
            IF (T.LT.0.0) THEN
               CALL LMINLS(X,W(IXIT+1),N,FMIN,M,FVECC,W(IXTEMP+1),
     +                     CALFUN,JFAIL)
               IF (JFAIL.EQ.1) IFAIL=1
               DO 240 J=1,N
                  W(IXI+N*(J-1)+IBIG)=W(IXIT+J)
240            CONTINUE
            END IF
         END IF
500   CONTINUE
      IFLAG=-2
      RETURN
      END