      SUBROUTINE DGSTPD(NR,N,G,A,P,SX,RNWTLN,DLT,NWTAKE,FSTDOG,
     +     SSD,V,CLN,ETA,SC,IPR,STEPMX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C PURPOSE
C -------
C FIND NEW STEP BY DOUBLE DOGLEG ALGORITHM
C
C PARAMETERS
C ----------
C NR           --> ROW DIMENSION OF MATRIX
C N            --> DIMENSION OF PROBLEM
C G(N)         --> GRADIENT AT CURRENT ITERATE, G(X)
C A(N,N)       --> CHOLESKY DECOMPOSITION OF HESSIAN IN
C                  LOWER PART AND DIAGONAL
C P(N)         --> NEWTON STEP
C SX(N)        --> DIAGONAL SCALING MATRIX FOR X
C RNWTLN       --> NEWTON STEP LENGTH
C DLT         <--> TRUST REGION RADIUS
C NWTAKE      <--> BOOLEAN, =.TRUE. IF NEWTON STEP TAKEN
C FSTDOG      <--> BOOLEAN, =.TRUE. IF ON FIRST LEG OF DOGLEG
C SSD(N)      <--> WORKSPACE [CAUCHY STEP TO THE MINIMUM OF THE
C                  QUADRATIC MODEL IN THE SCALED STEEPEST DESCENT
C                  DIRECTION] [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C V(N)        <--> WORKSPACE  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C CLN         <--> CAUCHY LENGTH
C                  [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C ETA              [RETAIN VALUE BETWEEN SUCCESSIVE CALLS]
C SC(N)       <--  CURRENT STEP
C IPR          --> DEVICE TO WHICH TO SEND OUTPUT
C STEPMX       --> MAXIMUM ALLOWABLE STEP SIZE
C
C INTERNAL VARIABLES
C ------------------
C CLN              LENGTH OF CAUCHY STEP
C
      DIMENSION G(N),P(N)
      DIMENSION SX(N)
      DIMENSION SC(N),SSD(N),V(N)
      DIMENSION A(NR,1)
      LOGICAL NWTAKE,FSTDOG
      IPR=IPR
C
C CAN WE TAKE NEWTON STEP
C
      IF(RNWTLN.GT.DLT) GO TO 100
C     IF(RNWTLN.LE.DLT)
C     THEN
        NWTAKE=.TRUE.
        DO 10 I=1,N
          SC(I)=P(I)
   10   CONTINUE
        DLT=RNWTLN
C$      WRITE(IPR,951)
        GO TO 700
C     ELSE
C
C NEWTON STEP TOO LONG
C CAUCHY STEP IS ON DOUBLE DOGLEG CURVE
C
  100   NWTAKE=.FALSE.
        IF(.NOT.FSTDOG) GO TO 200
C       IF(FSTDOG)
C       THEN
C
C         CALCULATE DOUBLE DOGLEG CURVE (SSD)
          FSTDOG=.FALSE.
          ALPHA=0.D0
          DO 110 I=1,N
            ALPHA=ALPHA + (G(I)*G(I))/(SX(I)*SX(I))
  110     CONTINUE
          BETA=0.D0
          DO 130 I=1,N
            TMP=0.D0
            DO 120 J=I,N
              TMP=TMP + (A(J,I)*G(J))/(SX(J)*SX(J))
  120       CONTINUE
            BETA=BETA+TMP*TMP
  130     CONTINUE
          DO 140 I=1,N
            SSD(I)=-(ALPHA/BETA)*G(I)/SX(I)
  140     CONTINUE
          CLN=ALPHA*SQRT(ALPHA)/BETA
          ETA=.2D0 + (.8D0*ALPHA*ALPHA)/(-BETA*DDOT(N,G,1,P,1))
          DO 150 I=1,N
            V(I)=ETA*SX(I)*P(I) - SSD(I)
  150     CONTINUE
          IF (DLT .EQ. (-1.0D0)) DLT = MIN(CLN, STEPMX)
C$        WRITE(IPR,954) ALPHA,BETA,CLN,ETA
C$        WRITE(IPR,955)
C$        WRITE(IPR,960) (SSD(I),I=1,N)
C$        WRITE(IPR,956)
C$        WRITE(IPR,960) (V(I),I=1,N)
C       ENDIF
  200   IF(ETA*RNWTLN.GT.DLT) GO TO 220
C       IF(ETA*RNWTLN .LE. DLT)
C       THEN
C
C         TAKE PARTIAL STEP IN NEWTON DIRECTION
C
          DO 210 I=1,N
            SC(I)=(DLT/RNWTLN)*P(I)
  210     CONTINUE
C$        WRITE(IPR,957)
          GO TO 700
C       ELSE
  220     IF(CLN.LT.DLT) GO TO 240
C         IF(CLN.GE.DLT)
C         THEN
C           TAKE STEP IN STEEPEST DESCENT DIRECTION
C
            DO 230 I=1,N
              SC(I)=(DLT/CLN)*SSD(I)/SX(I)
  230       CONTINUE
C$          WRITE(IPR,958)
            GO TO 700
C         ELSE
C
C           CALCULATE CONVEX COMBINATION OF SSD AND ETA*P
C           WHICH HAS SCALED LENGTH DLT
C
  240       DOT1=DDOT(N,V,1,SSD,1)
            DOT2=DDOT(N,V,1,V,1)
            ALAM=(-DOT1+SQRT((DOT1*DOT1)-DOT2*(CLN*CLN-DLT*DLT)))/DOT2
            DO 250 I=1,N
              SC(I)=(SSD(I) + ALAM*V(I))/SX(I)
  250       CONTINUE
C$          WRITE(IPR,959)
C         ENDIF
C       ENDIF
C     ENDIF
  700 CONTINUE
C$    WRITE(IPR,952) FSTDOG,NWTAKE,RNWTLN,DLT
C$    WRITE(IPR,953)
C$    WRITE(IPR,960) (SC(I),I=1,N)
      RETURN
C
  951 FORMAT(27H0DGSTPD    TAKE NEWTON STEP)
  952 FORMAT(18H DGSTPD    FSTDOG=,L1/
     +       18H DGSTPD    NWTAKE=,L1/
     +       18H DGSTPD    RNWTLN=,E20.13/
     +       18H DGSTPD    DLT   =,E20.13)
  953 FORMAT(28H DGSTPD    CURRENT STEP (SC))
  954 FORMAT(18H DGSTPD    ALPHA =,E20.13/
     +       18H DGSTPD    BETA  =,E20.13/
     +       18H DGSTPD    CLN   =,E20.13/
     +       18H DGSTPD    ETA   =,E20.13)
  955 FORMAT(28H DGSTPD    CAUCHY STEP (SSD))
  956 FORMAT(12H DGSTPD    V)
  957 FORMAT(48H0DGSTPD    TAKE PARTIAL STEP IN NEWTON DIRECTION)
  958 FORMAT(50H0DGSTPD    TAKE STEP IN STEEPEST DESCENT DIRECTION)
  959 FORMAT(39H0DGSTPD    TAKE CONVEX COMBINATION STEP)
  960 FORMAT(14H DGSTPD       ,5(E20.13,3X))
      END
