$NOFLOATCALLS
c$INCLUDE: 'DIST.H'
c$INCLUDE: '\FORTRAN\MATUTIL\MAT3.H'
      FUNCTION DIST(A,C)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 A(*),C(*), X
      INTEGER I
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FBASE, FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS,MOBSYM
      LOGICAL LNFLAG
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
      REAL*8 SCALE
      COMMON /SCALE/ SCALE(MAXELT)
      DIST=0.0
      DO 10 I=1,NELTS
         X=SCALE(I)*(A(I)-C(I))
         DIST=DIST+X*X
10    CONTINUE
      DIST=EXP(-DIST)
      RETURN
      END

      SUBROUTINE GRAD(A,GVEC)
C     Returns GVEC.  A is input, and is stored in A0.  A0 then
C     remains stable while various function evaluations are
C     made for a line search.
      IMPLICIT REAL*8(A-H,O-Z)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV,DET,FBASE,FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS,MOBSYM
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      LOGICAL LNFLAG
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
c----------------variables shared by GRAD and NEWT-----------
      INTEGER MELSYM
      PARAMETER(MELSYM=MAXELT*(MAXELT+1)/2)
      REAL*8 A0(MAXELT),G0VEC(MAXELT),H0(MELSYM),F0
      COMMON /TAYLOR/ A0,G0VEC,H0,F0
      SAVE /TAYLOR/
C------------------------------------------------------------------
      REAL*8 A(MAXELT),GVEC(MAXELT),CORVEC(MAXOBS),SCALE
     & ,SMTH, SSQ, DIST, DOT
      INTEGER I, IOBS
      COMMON /SCALE/ SCALE(MAXELT)
      CALL MATMOV(A,A0,NELTS)
      F0=SMTH(A0)
      DO 10 I=1,NELTS
         SSQ=SCALE(I)*SCALE(I)
         DO 11 IOBS=1,NOBS
11          CORVEC(IOBS) =-DIST(A,B(1,IOBS))*2.*SSQ
     1       *(A(I)-B(I,IOBS))
         GVEC(I)=DOT(CORVEC,TBFV,NOBS)
10    CONTINUE
      CALL MATMOV(GVEC,G0VEC,NELTS)
      RETURN
      END

      SUBROUTINE NEWT(A,ABORT)
      IMPLICIT REAL*8(A-H,O-Z)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FIT, FBASE
      INTEGER NELTS,NOBS,MAXELT,MAXOBS,MOBSYM
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      LOGICAL LNFLAG
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
      REAL*8 CORVEC(MAXELT),WORVEC(MAXELT),A(*), SCALE
     & ,FTOP, DOT, STEP0, STEP, F1, SMTH
      INTEGER IELT, J
      LOGICAL ANSYES,BOMB,ABORT
      COMMON /SCALE/ SCALE(MAXELT)
c----------------variables shared by GRAD and NEWT-----------
      INTEGER MELSYM
      PARAMETER(MELSYM=MAXELT*(MAXELT+1)/2)
      REAL*8 A0(MAXELT),G0VEC(MAXELT),H0(MELSYM),F0
      COMMON /TAYLOR/ A0,G0VEC,H0,F0
      SAVE /TAYLOR/
C-----------------------------------------------------------------
C           here WORVEC and CORVEC are being used for local storage,
C           without regard to their names.  WORVEC is minus the Gauss-Newton
C           step, and CORVEC is the trial coefficient vector.
C-----------------------------------------------
      CALL TRUNWT(WORVEC,BOMB)
      IF(BOMB) THEN
         WRITE(*,*)'Guess the function value achievable in this step: '
         READ (*,*) FTOP
         CALL MATMOV(G0VEC,WORVEC,NELTS)
C        The scaling below is to use diagonal of hessian as if there
C        were no off-diagonal elements.  This avoids dominance of the
C        gradient by b's with naturally small scales.
         DO 51 IELT=1,NELTS
51          WORVEC(IELT)=WORVEC(IELT)/(SCALE(IELT)*SCALE(IELT))
         CALL SCALAR(WORVEC,(FTOP-F0)/DOT(WORVEC,G0VEC,NELTS),NELTS)
      ENDIF
      STEP0=1.
      STEP=1.
50    CONTINUE
        CALL SCALAR(WORVEC,STEP/STEP0,NELTS)
        CALL MATADD(A0,WORVEC,CORVEC,NELTS)
        F1=SMTH(CORVEC)
        WRITE(*,*) 'Coeff. vector, function value: '
        IF(LNFLAG) THEN
            WRITE(*,*) (EXP(CORVEC(J)),J=1,NELTS), F1
        ELSE
            WRITE(*,*) (CORVEC(J),J=1,NELTS),F1
        ENDIF
        WRITE(*,*) 'Current step size is ', STEP,
     1 '.  Try new step size?'
        IF (ANSYES()) THEN
           STEP0=STEP
           READ(*,*) STEP
           GO TO 50
        ENDIF
      WRITE(*,*) 'Reset current coeff. vector? '
      IF (ANSYES()) THEN
        CALL MATADD(A0,WORVEC,A,NELTS)
        ABORT=.FALSE.
        RETURN
      ELSE
        ABORT=.TRUE.
        RETURN
      ENDIF
      END

      SUBROUTINE TRUNWT(DIREC,BOMB)
      IMPLICIT REAL*8(A-H,O-Z)
      LOGICAL BOMB
      REAL*8 DIREC(*)
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FBASE, FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS, MOBSYM
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      LOGICAL LNFLAG
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C---------------------------------------------
      REAL*8 CORVEC(MAXOBS),WORVEC(MAXOBS), SCALE, DIST, WTIJ, DOT
      INTEGER IOBS, II, IELT, JELT
      COMMON /SCALE/ SCALE(MAXELT)
c----------------variables shared by GRAD and NEWT-----------
      INTEGER MELSYM
      PARAMETER(MELSYM=MAXELT*(MAXELT+1)/2)
      REAL*8 A0(MAXELT),G0VEC(MAXELT),H0(MELSYM),F0
      COMMON /TAYLOR/ A0,G0VEC,H0,F0
         DO 13 IOBS=1,NOBS
13          WORVEC(IOBS)=DIST(B(1,IOBS),A0)
         II=1
         DO 10 IELT=1,NELTS
            DO 11 JELT=1,IELT
               WTIJ=(SCALE(IELT)*SCALE(JELT))
               WTIJ=WTIJ*WTIJ
               DO 12 IOBS=1,NOBS
                  CORVEC(IOBS)= -4.*(B(JELT,IOBS)-A0(JELT))*
     1              (B(IELT,IOBS)-A0(IELT))*WORVEC(IOBS)*WTIJ
12             CONTINUE
               IF(IELT.EQ.JELT) THEN
                  DO 14 IOBS=1,NOBS
                     CORVEC(IOBS)=CORVEC(IOBS)+
     1                  2.*SCALE(IELT)*SCALE(IELT)*WORVEC(IOBS)
14                CONTINUE
               ENDIF
               H0(II)=DOT(CORVEC,TBFV,NOBS)
               II=II+1
11          CONTINUE
10       CONTINUE
         CALL SYMINV(H0,NELTS,BOMB)
         IF (.NOT.BOMB) THEN
            CALL SYMVEC(H0,G0VEC,DIREC,NELTS)
         ELSE
            WRITE(*,*) 'Hessian not p.d.'
         ENDIF
         RETURN
         END

      SUBROUTINE AUTNWT(A,F,FINC,ABORT)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 A(*), FINC
      LOGICAL UP,DOWN,ABORT,BOMB, DONE
C     A returns with new parameter vector.
C--------------------------------------------
C     If a "natural increment" for the function is not around 20,
C     FINC, used to pick a step size in pure gradient search, should
C     be changed.
C----------------shared variables------------
      REAL*8 B,TBCVNV,FVAL,TBFV, DET, FBASE, FIT
      INTEGER NELTS,NOBS,MAXELT,MAXOBS, MOBSYM
      PARAMETER(MAXELT=15,MAXOBS=120,MOBSYM=MAXOBS*(MAXOBS+1)/2)
      LOGICAL LNFLAG
      COMMON /WORK/ DET, B(MAXELT,MAXOBS),TBCVNV(MOBSYM),
     1       FBASE,FVAL(MAXOBS),TBFV(MAXOBS),NELTS,NOBS,FIT,LNFLAG
C----------------variables shared by GRAD and NEWT-----------
      INTEGER MELSYM
      PARAMETER(MELSYM=MAXELT*(MAXELT+1)/2)
      REAL*8 A0(MAXELT),G0VEC(MAXELT),H0(MELSYM),F0
      COMMON /TAYLOR/ A0,G0VEC,H0,F0
C---------------------------------------------
      REAL*8 CORVEC(MAXELT),WORVEC(MAXELT),DOT,MINANG,PRECIS,SCALE
     & ,F1, F, SMTH, GSIZ, DDER, DDER0, WSIZ, STEP0, STEP, FACTOR
      INTEGER IELT
      PARAMETER (MINANG=.0001, PRECIS=1.D-30)
      COMMON /SCALE/ SCALE(MAXELT)
C-----------------------------------------------------------------
      ABORT=.FALSE.
      DONE = .FALSE.
      F1=SMTH(A)
      F=F1
      GSIZ = DOT(G0VEC,G0VEC,NELTS)
      IF (GSIZ .LE. PRECIS) THEN
         F=F1
         ABORT=.FALSE.
         WRITE(*,*) '0 gradient.'
         RETURN
      ENDIF
      CALL TRUNWT(WORVEC,BOMB)
      IF (BOMB) THEN
         CALL MATMOV(G0VEC,WORVEC,NELTS)
         DO 51 IELT=1,NELTS
51          WORVEC(IELT)=WORVEC(IELT)/(SCALE(IELT)*SCALE(IELT))
         DDER0=DOT(G0VEC,WORVEC,NELTS)
         CALL SCALAR(WORVEC,FINC/DDER0 , NELTS)
         DDER0=FINC
      ELSE
         DDER0 = DOT(WORVEC,G0VEC,NELTS)
         WSIZ=DOT(WORVEC,WORVEC,NELTS)
         IF( DDER0/SQRT(GSIZ*WSIZ) .LT. MINANG) THEN
            CALL MATMOV(G0VEC,CORVEC,NELTS)
            CALL SCALAR(CORVEC, .01*SQRT(WSIZ/GSIZ),NELTS)
            CALL MATADD(WORVEC,CORVEC,WORVEC,NELTS)
            DDER0=DOT(WORVEC,G0VEC,NELTS)
         ENDIF
      ENDIF
      if (dder0.lt.0) then
         write(*,*) 'bad gradient'
         stop
      endif
      STEP0=1.
      STEP=1.
      FACTOR = 3.
      UP = .TRUE.
      DOWN = .TRUE.
50    CONTINUE
        DDER = STEP*DDER0
        CALL SCALAR(WORVEC,STEP/STEP0,NELTS)
        CALL MATADD(A0,WORVEC,CORVEC,NELTS)
        F1=SMTH(CORVEC)
        IF (F1.GT.F) THEN 
            F=F1
            CALL MATMOV(CORVEC,A,NELTS)
        ENDIF
        IF ((F1-F0) .GT. .7*DDER) THEN
            IF (.NOT. UP) THEN
               FACTOR = FACTOR**.3
               IF (FACTOR .LT. 1.001) ABORT=.TRUE.
               UP=.TRUE.
            ENDIF
            DOWN=.FALSE.
            STEP0=STEP
            STEP = FACTOR*STEP
            IF(STEP.GT.1./PRECIS) ABORT=.TRUE.
        ELSE IF (F1-F0 .LT. .3*DDER) THEN
            IF(.NOT.DOWN) THEN
               FACTOR = FACTOR**.3
               IF( FACTOR .LT. 1.001) ABORT=.TRUE.
               DOWN=.TRUE.
            ENDIF
            UP = .FALSE.
            STEP0=STEP
            STEP = STEP/FACTOR
            IF(STEP.LT.PRECIS) ABORT=.TRUE.
        ELSE
            DONE = .TRUE.
        ENDIF
        IF (.NOT.(ABORT.OR.DONE)) THEN
            GO TO 50
        ENDIF
      FINC=F-F0
      RETURN
      END

      


