C
C*************************
C**       LSQMIN        **
C*************************
C
C
      SUBROUTINE LSQMIN(M,N,X,SHFLIM,MAXITR,H,DAMP,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          The routine calculates approximate partial derivatives using
C          the given step size H. A non-linear least squares refinement
C          is carried out.
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 SHFLIM (I)    If the absolute value of each shift is less
C                                than SHFT, then the refinement is assumed
C                                to be completed.
C INTEGER          MAXITR (I)    The maximum number of iterations to be taken
C DOUBLE PRECISION      H (I)    The step size to be used in calculating the
C                                partial derivatives.
C DOUBLE PRECISION   DAMP (I)    A damping factor to apply to the calculated
C                                shifts <= 1.0 (e.g. 0.75)
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 refinement not converged after MAXITR
C                                      iterations
C                                   -2 Normal matrix is singular
C DOUBLE PRECISION      W (*)    An array for work space dimensioned to at
C                                least (M+2*N+3)*N + 2*M
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 SHFLIM, FMIN, H, DAMP
      EXTERNAL CALFUN

      DOUBLE PRECISION W(N*N+4*N)
C
C Work array split up as follows:
C
C     XT(N)          Temporary copy of parameters
C     FC(M)          Function values for current values
C     FT(M)          Temp function values
C     A(M,N)         M*N derivs. matrix |A|
C     AN(N,N)        N*N Normal matrix
C     AI(N,N)        Inverse of the Normal matrix
C     RHS(N)         Right hand side elements
C     SHF(N)         Shifts
C
C
C====== LOCALS:
C
      LOGICAL CONVG
      DOUBLE PRECISION FSAV, DET
C
C-------------------------------------------------------------------------------
C
C
C====== Set offsets to the work array
C
      IXT=0
      IFC=IXT+N
      IFT=IFC+M
      IA=IFT+M
      IAN=IA+M*N
      IAI=IAN+N*N
      IRHS=IAI+N*N
      ISHF=IRHS+N
C
C====== Calculate the function value with the starting parameters
C
      CALL CALFUN(M,N,W(IFC+1),X)
      FMIN=0.0
      DO 20 I=1,M
         FMIN=FMIN+W(IFC+I)*W(IFC+I)
20    CONTINUE
      FSAV=FMIN
C
C====== Start Iterations Loop
C====== =====================
C
      DO 600 ITER=1,MAXITR
C
C====== Make a temporary copy of the parameters
C
         DO 10 I=1,N
            W(IXT+I)=X(I)
10       CONTINUE
C
C====== Calculate approximate partial derivatives and form |A| matrix
C
         II=IA
         DO 100 I=1,N
            W(IXT+I)=W(IXT+I)+H
            CALL CALFUN(M,N,W(IFT+1),W(IXT+1))
            DO 50 J=1,M
               II=II+1
               W(II)=(W(IFT+J)-W(IFC+J))/H
50          CONTINUE
            W(IXT+I)=X(I)
100      CONTINUE
C
C====== Form the Normal matrix of least squares |AT| * |A|
C
         IN=IAN
         DO 250 J=1,N
            I1=IA
            DO 200 I=1,N
               J1=IA+(J-1)*M
               IN=IN+1
               W(IN)=0.0
               DO 150 K=1,M
                  I1=I1+1
                  J1=J1+1
                  W(IN)=W(IN)+W(I1)*W(J1)
150            CONTINUE
200         CONTINUE
250      CONTINUE
C
C====== Form the right hand side elements |AT| * g
C
         IR=IRHS
         II=IA
         DO 350 I=1,N
            IR=IR+1
            W(IR)=0.0
            JJ=IFC
            DO 300 J=1,M
               II=II+1
               JJ=JJ+1
               W(IR)=W(IR)-W(II)*W(JJ)
300         CONTINUE
350         CONTINUE
C
C====== Invert the Normal matrix
C
         CALL INVMAT(W(IAN+1),W(IAI+1),N,N,0,DET,IERR)
         IF (IERR.EQ.1) THEN
            IFLAG=-2
            RETURN
         END IF
C
C====== Calculate the shifts
C
         IS=ISHF
         II=IAI
         DO 450 I=1,N
            IS=IS+1
            W(IS)=0.0
            JJ=IRHS
            DO 400 J=1,N
               II=II+1
               JJ=JJ+1
               W(IS)=W(IS)+W(II)*W(JJ)
400         CONTINUE
450      CONTINUE
C
C====== Apply the (damped) shifts
C
         IS=ISHF
         DO 500 I=1,N
            IS=IS+1
            X(I)=X(I)+DAMP*W(IS)
500      CONTINUE
C
C====== Calculate the function values with the shifted parameters
C
         CALL CALFUN(M,N,FVECC,X)
         FMIN=0.0
         DO 520 I=1,M
            FMIN=FMIN+FVECC(I)*FVECC(I)
            W(IFC+I)=FVECC(I)
520      CONTINUE
         IF (FMIN.GE.FSAV) THEN
            IFLAG=ITER
            FMIN=FSAV
            DO 525 I=1,N
               X(I)=W(IXT+I)
525         CONTINUE
            RETURN
         ENDIF
         CONVG=.TRUE.
         IS=ISHF
         DO 530 I=1,N
            IS=IS+1
            IF (ABS(W(IS)).GT.SHFLIM) THEN
               CONVG=.FALSE.
               GO TO 540
            END IF
530      CONTINUE
540      IF (CONVG) THEN
            IFLAG=ITER
            RETURN
         END IF
         FSAV=FMIN
600   CONTINUE
      IFLAG=-1
      RETURN
      END