C  Comments, bug reports, etc. should be sent to:  portnoy@stat.uiuc.edu
C
      SUBROUTINE CRQF(M,N,MPLUS,N2,X,Y,C,TZERO,MAXW,STEP,IFT,H,XH,
     * WA,WB,WC,WD,WE,WF,IA,NSOL,SOL,LSOL,ICEN,TCEN,LCEN)
C
C     M = number of Observations
C     N = Number of Parameters
C     MPLUS = row dim of X, WA, and WC .GE. M+1
C     N2 = N+2
C     X is the X matrix (MPLUS BY N) (includes censored obs)
C     Y is the y vector (M) (includes censored obs)
C     C is the censoring integer(M) vector, 1 if obs is censored,
C                                           0 otherwise
C     TZERO is the initial tau value
C     MAXW ( .LE. MPLUS ): Maximal number of calls to weighted RQ
C           for possible degeneracies. If exceeded (IFT = 7), dither.
C           If  MAXW .LE. 0, don't pivot -- use only weighted rq
C     STEP: Step size for weighted rq at degeneracy
C     IFT exit code:
C        0:  ok
C        1:  M < 0  or  N < 0  OR  M < N
C	 2:  MPLUS .LT. M + 1  or  N2 .NE. N+2
C        3:  MAXW > MPLUS
C        4:  less than N noncensored obs above the tau=0 soln
C        5:  possible degeneracy, rq called at tau + step, see IA
C        6:  LSOL becomes greater than NSOL
C        7:  MAXW exceeded
C        8:  weighted rq tries to include infinite basis element
C     H  is an integer work vector (N) ( = basis indices )
C     XH is a double precision work array (N by N)
C     WA is a double precision work array (MPLUS by N)
C     WB is a double precision work vector (MPLUS)
C     WC is a double precision work array (MPLUS by N+2)
C     WD is a double precision work vector (MPLUS)
C     WE is a double precision work vector (MPLUS)
C     WF is a double precision work vector (N)
C     IA is an integer work vector (MPLUS) 
C     NSOL is an (estimated) row dimension of the solution array
C          if NSOL < M, solutions returned only at tau = i/(NSOL-1)
C          if all solutions are desired, NSOL = 3*M is a good choice
C  on output:
C     SOL is the coefficient solution array (N+2 by NSOL)
C         first and last rows give tau bkpts and quantile
C         rows 2:(N+1) give the beta coefficients
C     LSOL is the actual dimension of the solution arrays
C          LSOL = NSOL  if NSOL < M, else LSOL .LE. NSOL
C     ICEN (M) indicates status of censored observations
C          = 0 if uncensored (or uncrossed censored while pivoting)
C          = 1 if crossed censored
C          = 2 if deleted as below tau = 0 solution
C          = 3 if above last (maximal tau) solution
C     TCEN (M) are the corresponding tau value where censored obs is
C          firstcrossed:  weight = (tau - TCEN(I))/(1 - TCEN(I))
C          TCEN = 1 if censored obs is never crossed (ICEN = 3)
C          TCEN = 0 if censored obs is deleted (ICEN = 2)
C     LCEN = number of censored observations
C     WD = list of first MPLUS tau values at which degeneracy
C          may have occurred (and weighted RQ was called)
C     H(1) = number of calls to weighted RQ (apparent degeneracy)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER H(N),ICEN(M),C(M),IA(MPLUS)
      LOGICAL APC, DEG
      DOUBLE PRECISION ONE
      DIMENSION X(M,N),Y(M),SOL(N2,NSOL),WC(MPLUS,N2),XH(N,N)
      DIMENSION WA(MPLUS,N),WB(MPLUS),WD(MPLUS),WE(MPLUS)
      DIMENSION F(2),WF(N),TCEN(M)
      DATA BIG/1.D17/
      DATA ZERO/0.00D0/
      DATA HALF/0.50D0/
      DATA ONE/1.00D0/
      DATA TWO/2.00D0/
C
C  CHECK DIMENSION PARAMETERS
C
      DEG = .FALSE.
      IFT=0
      N1 = NSOL-1
      LCEN = 0
      DO 2 I = 1,M
      IF(C(I) .EQ. 1) LCEN = LCEN + 1
   2  CONTINUE
      IF(MPLUS .LT. M+1) IFT = 2
      IF(N2 .NE. N+2) IFT = 2
      IF(M .LE. ZERO .OR. N. LE. ZERO .OR. M .LT. N) IFT = 1
      IF(MAXW .GT. MPLUS) IFT = 3
      IF(IFT .NE. 0) RETURN
C
C  INITIALIZATION
C
      TOLER = 1.D-11
      TOL1 = 10.0D0*TOLER
      IF(TZERO .LE. ZERO  .OR.  TZERO .GT. STEP) TZERO = STEP
      DO 5 I=1,M
        ICEN(I) = 0
        TCEN(I) = ONE
        WB(I) = Y(I)
        DO 5 J=1,N
   5    WA(I,J) = X(I,J)
      MA = M
C
C  GET TAU = O SOLUTION
C  DELETE CENSORED OBS BELOW TAU = 0 SOLN
C
  10  IF(MA .LE. N) THEN
        IFT = 4
        RETURN
      ENDIF
      CALL RQ1(MA,N,MPLUS,N2,WA,WB,TZERO,TOLER,IFT1,
     * WF,WE,IA,WC,WD)
      K = 0
      KL = 0
C
C  CHECK FOR CENSORED OBS BELOW TAU = 0 SOLN
C
      L = 0
      DO 20 I=1,M
      IF(ICEN(I) .EQ. 2) GO TO 20
      L = L + 1
      IF( C(I) .EQ. 1 .AND. WE(L) .LE. TOL1 ) THEN
        K = K + 1
        ICEN(I) = 2
        GO TO 20
      ENDIF
      IF( DABS(WE(L)) .GE. TOL1) GO TO 20
      KL = KL + 1
      IF(KL .LE. N) H(KL) = I
  20  CONTINUE
      IF( K .EQ. 0) GO TO 40
      MA = MA - K
      L = 0
C
C  DELETE CENSORED OBS BELOW TAU = 0 SOLN, AND RETRY
C 
      DO 30 I=1,M
        IF(ICEN(I) .EQ. 2) GO TO 30
        L = L + 1
        WB(L) = Y(I)
        DO 25 J=1,N
  25      WA(L,J) = X(I,J)
  30    CONTINUE
      GO TO 10
  40  CONTINUE
C
C  SAVE INITIAL SOLUTION, AND INITIALIZE PIVOTING LOOP
C
      SOL(1,1) = ZERO
      DO 50 J = 1,N
  50  SOL(J+1,1) = WF(J)
      NWRQ = 0
      APC = .FALSE.
      LSOL = 2
      TAU = TZERO
C
C  COMPUTE NEXT TAU
C  FIRST CHECK FOR DEGENERACY AT TZERO
C
      IF (KL .GT. N)  GO TO 500
C
C  GET X(H,)
C
      DO 70 K = 1,N
        I= H(K)
        DO 60 J = 1,N
  60    XH(K,J) = X(I,J)
  70  CONTINUE
C
C  GET X(H,) INVERSE
C
      CALL DGECO(XH,N,N,IA,V,WC)
      IF(V. LT. TOLER) GO TO 500
      JOB = 01
      CALL DGEDI(XH,N,N,IA,F,WC,JOB)
C
C  GET RESIDUALS
C
      DO 90 I=1,M
      S = Y(I)
      DO 80 J = 1,N
  80  S = S - WF(J)*X(I,J)
  90  WE(I) = S
C
C  BEGIN PIVOTING; WD = GRAD UPPER BD, IA = SIGN(GRAD DENOM)
C
 200  CONTINUE
      CALL GRAD(X,M,N,H,ICEN,TCEN,XH,WE,TOLER,IA,WC,WD)
      KL = 0
      KM = 1
      S = WD(1)
      IF(N .EQ. 1) GO TO 230
      DO 210 J=2,N
 210  S = DMIN1(WD(J),S)
      DO 220 J=1,N
      IF(WD(J) .GE. S + TOLER) GO TO 220
        KL = KL + 1
        KM = J
 220  CONTINUE
 230  CONTINUE
C     
C  IF ALL POS RESIDS CENSORED, RETURN; ELSE CHECK FOR DEGEN.
C     
      IF(APC) THEN
        SOL(1,LSOL) = DMAX1(S, TAU)
        DO 240 J=1,N
 240    SOL(J+1,LSOL) = WF(J)
        GO TO 600
      ENDIF
C
C  IF DEGENERACY, CALL WEIGHTED RQ
C
      IF (KL .GT. 1) GO TO 500
C     
C  CHECK FOR INFEASIBILITY CAUSED BY REWEIGHTING
C 
      IF(S .LT. TAU + TOL1) THEN
        LSOL = LSOL - 1
        TAU = TAU - .8*STEP
        GO TO 500
      ENDIF
      TAU = S
C
C  GET NEW OBSERVATION TO ENTER BASIS
C  FIRST DEFINE WD = BASIS INDICATOR = 1 IF I IN H(J)
C
      DO 250 I=1,M
 250  WD(I) = ZERO
      DO 260 J=1,N
 260  WD(H(J)) = ONE
      KN = 0
      D = BIG 
      KIN = 0
      DO 300 I=1,M
      IF(ICEN(I) .EQ. 2) GO TO 300
      S = ZERO
      DO 270 J=1,N
 270  S = S + X(I,J)*XH(J,KM) 
      S = IA(KM)*S
      IF(DABS(S).LT.TOL1 .OR. (C(I).EQ.1 .AND. ICEN(I).NE.1)
     *    .OR. WD(I).EQ.ONE) GO TO 300
      S = WE(I)/S
      IF(S .LT. TOL1) GO TO 300
      IF(S .LE. D - TOL1) THEN
        D = S
        KIN = I
        KN = 0
      ELSE
        IF(S .LE. D + TOL1) KN = 1
      ENDIF
 300  CONTINUE
      IF(KN .EQ. 1) GO TO 500
      H(KM) = KIN
C
C  UPDATE SOLUTION
C  GET NEW X(H,)
C
        DO 310 K = 1,N
        I = H(K)
        DO 310 J = 1,N
 310    XH(K,J) = X(I,J)
C
C  GET X(H,) INVERSE
C
 315  CALL DGECO(XH,N,N,IA,V,WA)
      IF(V. LT. TOLER) GO TO 500
      JOB = 01
      CALL DGEDI(XH,N,N,IA,F,WA,JOB)
C
C  GET BETA-HAT, RESIDS
C
      DO 340 K = 1,N
      S = ZERO
      DO 330 J = 1,N
 330  S = S + XH(K,J)*Y(H(J))
 340  WF(K) = S
 345  DO 360 I=1,M
      S = Y(I)
      DO 350 J = 1,N
 350  S = S - WF(J)*X(I,J)
 360  WE(I) = S
C
C  SAVE SOLUTION
C
 365  DO 380 I = (LSOL-1),N1
      IF(NSOL.GE.M .OR. TAU .GT. DFLOAT(I)/DFLOAT(N1) - 10*TOL1) THEN
        SOL(1,LSOL) = TAU
        DO 370 J = 1,N
 370    SOL(J+1,LSOL) = WF(J)
        LSOL = LSOL + 1
        GO TO 390
      ENDIF
 380  CONTINUE
 390  CONTINUE
C
C  CHECK FOR DIM(SOL) EXCEEDED
C
      IF(LSOL .GT. NSOL) THEN
        IFT = 6
        GO TO 600
      ENDIF
      IF(APC) THEN
        LSOL = LSOL - 1
        GO TO 600
      ENDIF
C
C  SET WT FOR CROSSED CENSORED OBSERVATIONS
C  APC = .TRUE.  IF ALL POS RESID ARE UNCROSSED CENSORED
C
      APC = .TRUE.
C  if the following are left uncommented: still get zero column
C      TAUW = TAU - STEP/TWO
C      IF(MAXW.GT.ZERO)THEN
C        TAUW=TAU
C      ENDIF
      DO 400 I=1,M
      IF(ICEN(I) .EQ. 2) GO TO 400
      IF(WE(I).GE.TOL1 .AND. (C(I).EQ.0 .OR.  ICEN(I).EQ.1))
     *   APC = .FALSE.
      IF(WE(I).LE.-TOL1 .AND. C(I).EQ.1 .AND. ICEN(I).EQ.0) THEN
C  
C at this point,  Y(I) is a crossed censored obs  (C_i)
C for the grid method  TAUW = TAU - STEP*(B2-Y(I))/(B2-B1)
C   where   B1 = x_i' beta_j   and   B2 = x_i' beta_(j+1)
C otherwise (for pivot), TAUW = TAU
C
        IF(MAXW.LT.0) THEN
         B1 = ZERO
         B2 = ZERO
         DO 396 J=1,N
          B1 = B1 + X(I,J) * SOL(J+1,LSOL - 2)
396       B2 = B2 + X(I,J) * SOL(J+1,LSOL - 1)
         A1 = (B2 - Y(I))/(B2-B1)
         TAUW = TAU - A1*STEP
        ELSE
         TAUW = TAU
        ENDIF
        ICEN(I) = 1
        TCEN(I) = TAUW
      ENDIF
 400  CONTINUE
      IF(APC) THEN
        IF(DEG) THEN
          IFT = 5
          LSOL = LSOL - 1
          GO TO 600
        ENDIF
        GO TO 200
      ENDIF
C
C  RETURN IF TAU > 1
C
      IF(TAU .GE. ONE - 10.0D0 * TOL1) GO TO 600
      IF(DEG) GO TO 500
      IF(MAXW .GT. 0) GO TO 200
C
C  POSSIBLE DEGENERACY -- USE WEIGHTED RQ1 TO RESOLVE
C
 500  NWRQ = NWRQ + 1
      IF(MAXW .GT. 0  .AND.  NWRQ .GT. MAXW) THEN
        IFT = 7
        LSOL = LSOL - 1
        GO TO 600
      ENDIF
      IF(NWRQ .LE. NSOL) SOL(N+2,NWRQ) = TAU
      TAU = TAU + STEP
      IF(TAU .GE. ONE) TAU = ONE - TOL1
      L = 0
      L1 = 0
      DO 506 J=1,N
 506  WF(J) = ZERO
      DO 530 I = 1,M
      IF(ICEN(I) .EQ. 0) THEN
        IF (C(I) .EQ. 0) THEN
          L1 = L1 + 1
          WB(L1) = Y(I)
          DO 510 J=1,N
 510      WA(L1,J) = X(I,J)
        ELSE
          DO 514 J = 1,N
 514      WF(J) = WF(J) + X(I,J)
        ENDIF
      ENDIF
      IF(ICEN(I) .EQ. 1) THEN
        L1 = L1 + 1
        L = L+1
        W = (TAU - TCEN(I))/(ONE - TCEN(I))
        WB(L1) = W * Y(I)
        DO 520 J = 1,N
          WA(L1,J) = W * X(I,J)
 520      WF(J) = WF(J) + (ONE - W) * X(I,J)
      ENDIF
 530  CONTINUE
      MAL = L1+1
      DO 534 J=1,N
 534  WA(MAL,J) = WF(J)
      WB(MAL) = BIG
      CALL RQ1(MAL,N,MPLUS,N2,WA,WB,TAU,TOLER,IFT1,WF,WE,IA,WC,WD)
      DEG = .FALSE.
      IF(DABS(WE(MAL)) .LE. .0001) THEN
        IFT = 8
        LSOL = LSOL - 1
        GO TO 600
      ENDIF 
      L = 0
      K = 0
      DO 550 I=1,M
        IF(ICEN(I) .EQ. 2) GO TO 550
        L = L+1
        IF(DABS(WE(L)) .GE. TOL1) GO TO 550
          K = K+1
          IF(K .LE. N) THEN
            H(K) = I
            DO 540 J=1,N
 540        XH(K,J) = X(I,J)
          ENDIF
 550  CONTINUE
      IF(K .LT. N) THEN
        IFT = 8
        LSOL = LSOL - 1
        GO TO 600
      ENDIF
      IF(K .GT. N) DEG = .TRUE.
      GO TO 345
C
C  DEFINE OUTPUT AND RETURN
C
 600  H(1) = NWRQ
      L = MIN0(NWRQ,MPLUS)
      DO 610 I=1,L
 610  WD(I) = SOL(N+2,I)
      DO 620 I=1,M
      IF(ICEN(I) .EQ. 2) GO TO 620
      IF(C(I) .EQ.1 .AND. TCEN(I) .EQ. ONE) ICEN(I) = 3
 620  CONTINUE
      V = DFLOAT(MA)
      DO 650 J = 1,N
      S = ZERO
        DO 640 I = 1,M
          IF(ICEN(I) .EQ. 2) GO TO 640
          S = S + X(I,J)
 640    CONTINUE
 650  WF(J) = S/V
      DO 670 I = 1,LSOL
      S = ZERO
        DO 660 J = 1,N
 660    S = S + WF(J) * SOL(J+1,I)
 670  SOL(N+2,I) = S
      RETURN
      END
      SUBROUTINE GRAD(X,M,N,H,ICEN,TCEN,XH,
     * R,TOL,IFLAG,WA,GUP)
C
C  X matrix (M BY N)
C  M = Number of Observations
C  N = Number of Parameters
C  H = basis, integer(N) vector 
C  ICEN (integer(M)) = 2 for deleted obs
C                    = 1 for crossed non-censored obs
C                    = 0 otherwise
C  TCEN(M) are the corresponding tau values where a censored obs is
C    first crossed (ICEN=1): weight = (tau - TCEN(I))/(1 - TCEN(I))
C    TCEN = 1 if censored obs is never crossed
C  XH = (N by N) X(H,) inverse matrix
C  R(M) = residuals
C  TOL = zero tolerance (1.D-10 by default)
C  IFLAG (M)  work vector
C  WA (M by N) work array
C  returns:
C    GUP(N)  = upper bounds on tau
C    IFLAG(1:N) = sign(grad denom)
C               = +1 if bound from lower grad bound
C               = -1 if bound from upper grad bound
C    
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER H(N),ICEN(M),IFLAG(*)
      DOUBLE PRECISION ONE
      DIMENSION X(M,N),GUP(N),XH(N,N)
      DIMENSION R(M),TCEN(M),WA(M,N)
      DATA  ZERO/0.00D0/, ONE/1.0d0/
C
C GET X'(XH-INVERSE)
C
      DO 80 I = 1,M
      IF(ICEN(I) .EQ. 2) GO TO 80
      DO 70 J = 1,N
      A = ZERO
        DO 60 K = 1,N
  60    A = A + X(I,K)*XH(K,J)
  70  WA(I,J) = A
  80  CONTINUE
  
C
C  GET GRADIENT BOUNDS
C  FIRST SET IFLAG = 1 FOR BASIS INDICES (TEMPORARY)
C
      T = ZERO
      DO 90 I = 1,M
  90  IFLAG(I) = 0
      DO 95 J = 1,N
  95  IFLAG(H(J)) = 1
      DO 120 J=1,N
      A = ZERO
      B = ZERO
      C = ZERO
      D = ZERO
      DO 100 I = 1,M
        IF(ICEN(I) .EQ. 2) GO TO 100
        IF(ICEN(I) .EQ. 0) THEN
          IF(R(I) .GT. TOL)    A = A + WA(I,J)
          IF(R(I) .LT. - TOL)  B = B + WA(I,J)
          GO TO 100
        ENDIF
        IF(IFLAG(I).EQ.1) GO TO 100
C
C  IFLAG = 0: NOT BASIS  AND  ICEN = 1: CROSSED CENSORED
C
        IF(R(I) .LT. -TOL) THEN
           T = TCEN(I)/(ONE - TCEN(I))
           C = C - T*WA(I,J)
           GO TO 100
        ENDIF
        IF(R(I) .GT. TOL) THEN          
           D = D - WA(I,J)
        ENDIF
 100  CONTINUE
      D = D - C
      L = ICEN(H(J))
      IF(L .NE. 0) T = TCEN(H(J))/( ONE - TCEN(H(J)) )
      S = DFLOAT(L)*(T + ONE) - ONE
      E1 = A + B - D - S
      E2 = A + B - D + ONE
      IF(E1 .GT. ZERO) THEN
        GUP(J) = (B + C - S)/E1
        IFLAG(J+M) = 1
        GO TO 120
      ENDIF
      IF(E2 .LT. ZERO) THEN
        GUP(J) = (B + C)/E2
        IFLAG(J+M) = -1
        GO TO 120
      ENDIF
      GUP(J) = - ONE
 120  CONTINUE
      DO 130 J = 1,N
 130  IFLAG(J) = IFLAG(J+M)
      RETURN
      END