1C  Comments, bug reports, etc. should be sent to:  portnoy@stat.uiuc.edu
2C
3      SUBROUTINE CRQF(M,N,MPLUS,N2,X,Y,C,TZERO,MAXW,STEP,IFT,H,XH,
4     * WA,WB,WC,WD,WE,WF,IA,NSOL,SOL,LSOL,ICEN,TCEN,LCEN)
5C
6C     M = number of Observations
7C     N = Number of Parameters
8C     MPLUS = row dim of X, WA, and WC .GE. M+1
9C     N2 = N+2
10C     X is the X matrix (MPLUS BY N) (includes censored obs)
11C     Y is the y vector (M) (includes censored obs)
12C     C is the censoring integer(M) vector, 1 if obs is censored,
13C                                           0 otherwise
14C     TZERO is the initial tau value
15C     MAXW ( .LE. MPLUS ): Maximal number of calls to weighted RQ
16C           for possible degeneracies. If exceeded (IFT = 7), dither.
17C           If  MAXW .LE. 0, don't pivot -- use only weighted rq
18C     STEP: Step size for weighted rq at degeneracy
19C     IFT exit code:
20C        0:  ok
21C        1:  M < 0  or  N < 0  OR  M < N
22C	 2:  MPLUS .LT. M + 1  or  N2 .NE. N+2
23C        3:  MAXW > MPLUS
24C        4:  less than N noncensored obs above the tau=0 soln
25C        5:  possible degeneracy, rq called at tau + step, see IA
26C        6:  LSOL becomes greater than NSOL
27C        7:  MAXW exceeded
28C        8:  weighted rq tries to include infinite basis element
29C     H  is an integer work vector (N) ( = basis indices )
30C     XH is a double precision work array (N by N)
31C     WA is a double precision work array (MPLUS by N)
32C     WB is a double precision work vector (MPLUS)
33C     WC is a double precision work array (MPLUS by N+2)
34C     WD is a double precision work vector (MPLUS)
35C     WE is a double precision work vector (MPLUS)
36C     WF is a double precision work vector (N)
37C     IA is an integer work vector (MPLUS)
38C     NSOL is an (estimated) row dimension of the solution array
39C          if NSOL < M, solutions returned only at tau = i/(NSOL-1)
40C          if all solutions are desired, NSOL = 3*M is a good choice
41C  on output:
42C     SOL is the coefficient solution array (N+2 by NSOL)
43C         first and last rows give tau bkpts and quantile
44C         rows 2:(N+1) give the beta coefficients
45C     LSOL is the actual dimension of the solution arrays
46C          LSOL = NSOL  if NSOL < M, else LSOL .LE. NSOL
47C     ICEN (M) indicates status of censored observations
48C          = 0 if uncensored (or uncrossed censored while pivoting)
49C          = 1 if crossed censored
50C          = 2 if deleted as below tau = 0 solution
51C          = 3 if above last (maximal tau) solution
52C     TCEN (M) are the corresponding tau value where censored obs is
53C          firstcrossed:  weight = (tau - TCEN(I))/(1 - TCEN(I))
54C          TCEN = 1 if censored obs is never crossed (ICEN = 3)
55C          TCEN = 0 if censored obs is deleted (ICEN = 2)
56C     LCEN = number of censored observations
57C     WD = list of first MPLUS tau values at which degeneracy
58C          may have occurred (and weighted RQ was called)
59C     H(1) = number of calls to weighted RQ (apparent degeneracy)
60C
61      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
62      INTEGER H(N),ICEN(M),C(M),IA(MPLUS)
63      LOGICAL APC, DEG
64      DOUBLE PRECISION ONE
65      DIMENSION X(M,N),Y(M),SOL(N2,NSOL),WC(MPLUS,N2),XH(N,N)
66      DIMENSION WA(MPLUS,N),WB(MPLUS),WD(MPLUS),WE(MPLUS)
67      DIMENSION F(2),WF(N),TCEN(M)
68      DATA BIG/1.D17/
69      DATA ZERO/0.00D0/
70      DATA HALF/0.50D0/
71      DATA ONE/1.00D0/
72      DATA TWO/2.00D0/
73C
74C  CHECK DIMENSION PARAMETERS
75C
76      DEG = .FALSE.
77      IFT=0
78      N1 = NSOL-1
79      LCEN = 0
80      DO 2 I = 1,M
81      IF(C(I) .EQ. 1) LCEN = LCEN + 1
82   2  CONTINUE
83      IF(MPLUS .LT. M+1) IFT = 2
84      IF(N2 .NE. N+2) IFT = 2
85      IF(M .LE. ZERO .OR. N. LE. ZERO .OR. M .LT. N) IFT = 1
86      IF(MAXW .GT. MPLUS) IFT = 3
87      IF(IFT .NE. 0) RETURN
88C
89C  INITIALIZATION
90C
91      TOLER = 1.D-11
92      TOL1 = 10.0D0*TOLER
93      IF(TZERO .LE. ZERO  .OR.  TZERO .GT. STEP) TZERO = STEP
94      DO 5 I=1,M
95        ICEN(I) = 0
96        TCEN(I) = ONE
97        WB(I) = Y(I)
98        DO 5 J=1,N
99   5    WA(I,J) = X(I,J)
100      MA = M
101C
102C  GET TAU = O SOLUTION
103C  DELETE CENSORED OBS BELOW TAU = 0 SOLN
104C
105  10  IF(MA .LE. N) THEN
106        IFT = 4
107        RETURN
108      ENDIF
109      CALL RQ1(MA,N,MPLUS,N2,WA,WB,TZERO,TOLER,IFT1,
110     * WF,WE,IA,WC,WD)
111      K = 0
112      KL = 0
113C
114C  CHECK FOR CENSORED OBS BELOW TAU = 0 SOLN
115C
116      L = 0
117      DO 20 I=1,M
118      IF(ICEN(I) .EQ. 2) GO TO 20
119      L = L + 1
120      IF( C(I) .EQ. 1 .AND. WE(L) .LE. TOL1 ) THEN
121        K = K + 1
122        ICEN(I) = 2
123        GO TO 20
124      ENDIF
125      IF( DABS(WE(L)) .GE. TOL1) GO TO 20
126      KL = KL + 1
127      IF(KL .LE. N) H(KL) = I
128  20  CONTINUE
129      IF( K .EQ. 0) GO TO 40
130      MA = MA - K
131      L = 0
132C
133C  DELETE CENSORED OBS BELOW TAU = 0 SOLN, AND RETRY
134C
135      DO 30 I=1,M
136        IF(ICEN(I) .EQ. 2) GO TO 30
137        L = L + 1
138        WB(L) = Y(I)
139        DO 25 J=1,N
140  25      WA(L,J) = X(I,J)
141  30    CONTINUE
142      GO TO 10
143  40  CONTINUE
144C
145C  SAVE INITIAL SOLUTION, AND INITIALIZE PIVOTING LOOP
146C
147      SOL(1,1) = ZERO
148      DO 50 J = 1,N
149  50  SOL(J+1,1) = WF(J)
150      NWRQ = 0
151      APC = .FALSE.
152      LSOL = 2
153      TAU = TZERO
154C
155C  COMPUTE NEXT TAU
156C  FIRST CHECK FOR DEGENERACY AT TZERO
157C
158      IF (KL .GT. N)  GO TO 500
159C
160C  GET X(H,)
161C
162      DO 70 K = 1,N
163        I= H(K)
164        DO 60 J = 1,N
165  60    XH(K,J) = X(I,J)
166  70  CONTINUE
167C
168C  GET X(H,) INVERSE
169C
170      CALL DGECO(XH,N,N,IA,V,WC)
171      IF(V. LT. TOLER) GO TO 500
172      JOB = 01
173      CALL DGEDI(XH,N,N,IA,F,WC,JOB)
174C
175C  GET RESIDUALS
176C
177      DO 90 I=1,M
178      S = Y(I)
179      DO 80 J = 1,N
180  80  S = S - WF(J)*X(I,J)
181  90  WE(I) = S
182C
183C  BEGIN PIVOTING; WD = GRAD UPPER BD, IA = SIGN(GRAD DENOM)
184C
185 200  CONTINUE
186      CALL GRAD(X,M,N,H,ICEN,TCEN,XH,WE,TOLER,IA,WC,WD)
187      KL = 0
188      KM = 1
189      S = WD(1)
190      IF(N .EQ. 1) GO TO 230
191      DO 210 J=2,N
192 210  S = DMIN1(WD(J),S)
193      DO 220 J=1,N
194      IF(WD(J) .GE. S + TOLER) GO TO 220
195        KL = KL + 1
196        KM = J
197 220  CONTINUE
198 230  CONTINUE
199C
200C  IF ALL POS RESIDS CENSORED, RETURN; ELSE CHECK FOR DEGEN.
201C
202      IF(APC) THEN
203        SOL(1,LSOL) = DMAX1(S, TAU)
204        DO 240 J=1,N
205 240    SOL(J+1,LSOL) = WF(J)
206        GO TO 600
207      ENDIF
208C
209C  IF DEGENERACY, CALL WEIGHTED RQ
210C
211      IF (KL .GT. 1) GO TO 500
212C
213C  CHECK FOR INFEASIBILITY CAUSED BY REWEIGHTING
214C
215      IF(S .LT. TAU + TOL1) THEN
216        LSOL = LSOL - 1
217        TAU = TAU - .8*STEP
218        GO TO 500
219      ENDIF
220      TAU = S
221C
222C  GET NEW OBSERVATION TO ENTER BASIS
223C  FIRST DEFINE WD = BASIS INDICATOR = 1 IF I IN H(J)
224C
225      DO 250 I=1,M
226 250  WD(I) = ZERO
227      DO 260 J=1,N
228 260  WD(H(J)) = ONE
229      KN = 0
230      D = BIG
231      KIN = 0
232      DO 300 I=1,M
233      IF(ICEN(I) .EQ. 2) GO TO 300
234      S = ZERO
235      DO 270 J=1,N
236 270  S = S + X(I,J)*XH(J,KM)
237      S = IA(KM)*S
238      IF(DABS(S).LT.TOL1 .OR. (C(I).EQ.1 .AND. ICEN(I).NE.1)
239     *    .OR. WD(I).EQ.ONE) GO TO 300
240      S = WE(I)/S
241      IF(S .LT. TOL1) GO TO 300
242      IF(S .LE. D - TOL1) THEN
243        D = S
244        KIN = I
245        KN = 0
246      ELSE
247        IF(S .LE. D + TOL1) KN = 1
248      ENDIF
249 300  CONTINUE
250      IF(KN .EQ. 1) GO TO 500
251      H(KM) = KIN
252C
253C  UPDATE SOLUTION
254C  GET NEW X(H,)
255C
256        DO 310 K = 1,N
257        I = H(K)
258        DO 310 J = 1,N
259 310    XH(K,J) = X(I,J)
260C
261C  GET X(H,) INVERSE
262C
263 315  CALL DGECO(XH,N,N,IA,V,WA)
264      IF(V. LT. TOLER) GO TO 500
265      JOB = 01
266      CALL DGEDI(XH,N,N,IA,F,WA,JOB)
267C
268C  GET BETA-HAT, RESIDS
269C
270      DO 340 K = 1,N
271      S = ZERO
272      DO 330 J = 1,N
273 330  S = S + XH(K,J)*Y(H(J))
274 340  WF(K) = S
275 345  DO 360 I=1,M
276      S = Y(I)
277      DO 350 J = 1,N
278 350  S = S - WF(J)*X(I,J)
279 360  WE(I) = S
280C
281C  SAVE SOLUTION
282C
283 365  DO 380 I = (LSOL-1),N1
284      IF(NSOL.GE.M .OR. TAU .GT. DFLOAT(I)/DFLOAT(N1) - 10*TOL1) THEN
285        SOL(1,LSOL) = TAU
286        DO 370 J = 1,N
287 370    SOL(J+1,LSOL) = WF(J)
288        LSOL = LSOL + 1
289        GO TO 390
290      ENDIF
291 380  CONTINUE
292 390  CONTINUE
293C
294C  CHECK FOR DIM(SOL) EXCEEDED
295C
296      IF(LSOL .GT. NSOL) THEN
297        IFT = 6
298        GO TO 600
299      ENDIF
300      IF(APC) THEN
301        LSOL = LSOL - 1
302        GO TO 600
303      ENDIF
304C
305C  SET WT FOR CROSSED CENSORED OBSERVATIONS
306C  APC = .TRUE.  IF ALL POS RESID ARE UNCROSSED CENSORED
307C
308      APC = .TRUE.
309C  if the following are left uncommented: still get zero column
310C      TAUW = TAU - STEP/TWO
311C      IF(MAXW.GT.ZERO)THEN
312C        TAUW=TAU
313C      ENDIF
314      DO 400 I=1,M
315      IF(ICEN(I) .EQ. 2) GO TO 400
316      IF(WE(I).GE.TOL1 .AND. (C(I).EQ.0 .OR.  ICEN(I).EQ.1))
317     *   APC = .FALSE.
318      IF(WE(I).LE.-TOL1 .AND. C(I).EQ.1 .AND. ICEN(I).EQ.0) THEN
319C
320C at this point,  Y(I) is a crossed censored obs  (C_i)
321C for the grid method  TAUW = TAU - STEP*(B2-Y(I))/(B2-B1)
322C   where   B1 = x_i' beta_j   and   B2 = x_i' beta_(j+1)
323C otherwise (for pivot), TAUW = TAU
324C
325        IF(MAXW.LT.0) THEN
326         B1 = ZERO
327         B2 = ZERO
328         DO 396 J=1,N
329          B1 = B1 + X(I,J) * SOL(J+1,LSOL - 2)
330396       B2 = B2 + X(I,J) * SOL(J+1,LSOL - 1)
331         A1 = (B2 - Y(I))/(B2-B1)
332         TAUW = TAU - A1*STEP
333        ELSE
334         TAUW = TAU
335        ENDIF
336        ICEN(I) = 1
337        TCEN(I) = TAUW
338      ENDIF
339 400  CONTINUE
340      IF(APC) THEN
341        IF(DEG) THEN
342          IFT = 5
343          LSOL = LSOL - 1
344          GO TO 600
345        ENDIF
346        GO TO 200
347      ENDIF
348C
349C  RETURN IF TAU > 1
350C
351      IF(TAU .GE. ONE - 10.0D0 * TOL1) GO TO 600
352      IF(DEG) GO TO 500
353      IF(MAXW .GT. 0) GO TO 200
354C
355C  POSSIBLE DEGENERACY -- USE WEIGHTED RQ1 TO RESOLVE
356C
357 500  NWRQ = NWRQ + 1
358      IF(MAXW .GT. 0  .AND.  NWRQ .GT. MAXW) THEN
359        IFT = 7
360        LSOL = LSOL - 1
361        GO TO 600
362      ENDIF
363      IF(NWRQ .LE. NSOL) SOL(N+2,NWRQ) = TAU
364      TAU = TAU + STEP
365      IF(TAU .GE. ONE) TAU = ONE - TOL1
366      L = 0
367      L1 = 0
368      DO 506 J=1,N
369 506  WF(J) = ZERO
370      DO 530 I = 1,M
371      IF(ICEN(I) .EQ. 0) THEN
372        IF (C(I) .EQ. 0) THEN
373          L1 = L1 + 1
374          WB(L1) = Y(I)
375          DO 510 J=1,N
376 510      WA(L1,J) = X(I,J)
377        ELSE
378          DO 514 J = 1,N
379 514      WF(J) = WF(J) + X(I,J)
380        ENDIF
381      ENDIF
382      IF(ICEN(I) .EQ. 1) THEN
383        L1 = L1 + 1
384        L = L+1
385        W = (TAU - TCEN(I))/(ONE - TCEN(I))
386        WB(L1) = W * Y(I)
387        DO 520 J = 1,N
388          WA(L1,J) = W * X(I,J)
389 520      WF(J) = WF(J) + (ONE - W) * X(I,J)
390      ENDIF
391 530  CONTINUE
392      MAL = L1+1
393      DO 534 J=1,N
394 534  WA(MAL,J) = WF(J)
395      WB(MAL) = BIG
396      CALL RQ1(MAL,N,MPLUS,N2,WA,WB,TAU,TOLER,IFT1,WF,WE,IA,WC,WD)
397      DEG = .FALSE.
398      IF(DABS(WE(MAL)) .LE. .0001) THEN
399        IFT = 8
400        LSOL = LSOL - 1
401        GO TO 600
402      ENDIF
403      L = 0
404      K = 0
405      DO 550 I=1,M
406        IF(ICEN(I) .EQ. 2) GO TO 550
407        L = L+1
408        IF(DABS(WE(L)) .GE. TOL1) GO TO 550
409          K = K+1
410          IF(K .LE. N) THEN
411            H(K) = I
412            DO 540 J=1,N
413 540        XH(K,J) = X(I,J)
414          ENDIF
415 550  CONTINUE
416      IF(K .LT. N) THEN
417        IFT = 8
418        LSOL = LSOL - 1
419        GO TO 600
420      ENDIF
421      IF(K .GT. N) DEG = .TRUE.
422      GO TO 345
423C
424C  DEFINE OUTPUT AND RETURN
425C
426 600  H(1) = NWRQ
427      L = MIN0(NWRQ,MPLUS)
428      DO 610 I=1,L
429 610  WD(I) = SOL(N+2,I)
430      DO 620 I=1,M
431      IF(ICEN(I) .EQ. 2) GO TO 620
432      IF(C(I) .EQ.1 .AND. TCEN(I) .EQ. ONE) ICEN(I) = 3
433 620  CONTINUE
434      V = DFLOAT(MA)
435      DO 650 J = 1,N
436      S = ZERO
437        DO 640 I = 1,M
438          IF(ICEN(I) .EQ. 2) GO TO 640
439          S = S + X(I,J)
440 640    CONTINUE
441 650  WF(J) = S/V
442      DO 670 I = 1,LSOL
443      S = ZERO
444        DO 660 J = 1,N
445 660    S = S + WF(J) * SOL(J+1,I)
446 670  SOL(N+2,I) = S
447      RETURN
448      END
449      SUBROUTINE GRAD(X,M,N,H,ICEN,TCEN,XH,
450     * R,TOL,IFLAG,WA,GUP)
451C
452C  X matrix (M BY N)
453C  M = Number of Observations
454C  N = Number of Parameters
455C  H = basis, integer(N) vector
456C  ICEN (integer(M)) = 2 for deleted obs
457C                    = 1 for crossed non-censored obs
458C                    = 0 otherwise
459C  TCEN(M) are the corresponding tau values where a censored obs is
460C    first crossed (ICEN=1): weight = (tau - TCEN(I))/(1 - TCEN(I))
461C    TCEN = 1 if censored obs is never crossed
462C  XH = (N by N) X(H,) inverse matrix
463C  R(M) = residuals
464C  TOL = zero tolerance (1.D-10 by default)
465C  IFLAG (M)  work vector
466C  WA (M by N) work array
467C  returns:
468C    GUP(N)  = upper bounds on tau
469C    IFLAG(1:N) = sign(grad denom)
470C               = +1 if bound from lower grad bound
471C               = -1 if bound from upper grad bound
472C
473      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
474      INTEGER H(N),ICEN(M),IFLAG(*)
475      DOUBLE PRECISION ONE
476      DIMENSION X(M,N),GUP(N),XH(N,N)
477      DIMENSION R(M),TCEN(M),WA(M,N)
478      DATA  ZERO/0.00D0/, ONE/1.0d0/
479C
480C GET X'(XH-INVERSE)
481C
482      DO 80 I = 1,M
483      IF(ICEN(I) .EQ. 2) GO TO 80
484      DO 70 J = 1,N
485      A = ZERO
486        DO 60 K = 1,N
487  60    A = A + X(I,K)*XH(K,J)
488  70  WA(I,J) = A
489  80  CONTINUE
490
491C
492C  GET GRADIENT BOUNDS
493C  FIRST SET IFLAG = 1 FOR BASIS INDICES (TEMPORARY)
494C
495      T = ZERO
496      DO 90 I = 1,M
497  90  IFLAG(I) = 0
498      DO 95 J = 1,N
499  95  IFLAG(H(J)) = 1
500      DO 120 J=1,N
501      A = ZERO
502      B = ZERO
503      C = ZERO
504      D = ZERO
505      DO 100 I = 1,M
506        IF(ICEN(I) .EQ. 2) GO TO 100
507        IF(ICEN(I) .EQ. 0) THEN
508          IF(R(I) .GT. TOL)    A = A + WA(I,J)
509          IF(R(I) .LT. - TOL)  B = B + WA(I,J)
510          GO TO 100
511        ENDIF
512        IF(IFLAG(I).EQ.1) GO TO 100
513C
514C  IFLAG = 0: NOT BASIS  AND  ICEN = 1: CROSSED CENSORED
515C
516        IF(R(I) .LT. -TOL) THEN
517           T = TCEN(I)/(ONE - TCEN(I))
518           C = C - T*WA(I,J)
519           GO TO 100
520        ENDIF
521        IF(R(I) .GT. TOL) THEN
522           D = D - WA(I,J)
523        ENDIF
524 100  CONTINUE
525      D = D - C
526      L = ICEN(H(J))
527      IF(L .NE. 0) T = TCEN(H(J))/( ONE - TCEN(H(J)) )
528      S = DFLOAT(L)*(T + ONE) - ONE
529      E1 = A + B - D - S
530      E2 = A + B - D + ONE
531      IF(E1 .GT. ZERO) THEN
532        GUP(J) = (B + C - S)/E1
533        IFLAG(J+M) = 1
534        GO TO 120
535      ENDIF
536      IF(E2 .LT. ZERO) THEN
537        GUP(J) = (B + C)/E2
538        IFLAG(J+M) = -1
539        GO TO 120
540      ENDIF
541      GUP(J) = - ONE
542 120  CONTINUE
543      DO 130 J = 1,N
544 130  IFLAG(J) = IFLAG(J+M)
545      RETURN
546      END
547