#     this software has been submitted to statlib and may be freely
#     used and redistributed for non-commercial purposes.  no guarantees
#     are offered or implied.  comments, bug reports, etc are welcome
#     and should be sent to roger@ysidro.econ.uiuc.edu or to
#
#				roger koenker
#				department of economics
#				university of illinois
#				champaign, illinois, 61820
#
#
subroutine rq(m,nn,m5,n2,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,ndsol,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1)
#
#     number of observations
#     number of parameters
#     m+5
#     n+2
#     n+4
#     a is the x matrix
#     b is the y vector
#     t, the desired quantile
#     toler, smallest detectable |x-y|/x machine precision to the 2/3
#     ift exit code:
#		0-ok
#		else dimensions inconsistent
#     x the parameter estimate betahat
#     e is the residual vector
#     s is an integer work array
#     wa is a double precision work array
#     wb is another double precision work array
#     nsol is an estimated (row) dimension of the primal solution array say 3*m; minimum value = 2
#     ndsol is an estimated (row) dimension of the dual solution array say 3*m; minimum value = 2
#     sol is the primal solution array
#     dsol is the dual solution arry
#     lsol is the actual dimension of the solution arrays
#     h is the matrix of basic observations indices
#     qn is the vector of residuals variances from the projection of each column of the x matrix on the remaining columns
#     cutoff, the critical point for N(0,1)
#     ci is the matrix of confident intervals
#     tnmat is the matrix of the JGPK rank test statistics
#     big, large positive finite floating-point number
#     utilization:  if you just want a solution at a single quantile you
#     needn't bother with sol, nsol, etc, if you want all the solutions
#     then set theta to something <0 and sol and dsol will return all the
#     estimated quantile solutions.
#     the algorithm  is a slightly modified version of algorithm as 229
#     described in koenker and d'orey, "computing regression quantiles,
#     applied statistics, pp. 383-393.
#
integer i,j,k,kl,kount,kr,l,lsol,m,m1,m2,m3,m4,m5,ift
integer n,n1,n2,nsol,ndsol,out,s(m),h(nn,nsol)
integer nn,n4,idxcf
logical stage,test,init,iend,lup
logical lci1,lci2,skip
double precision a1,aux,b1,big,d,dif,pivot,smax,t,t0,t1,tnt
double precision min,max,toler,zero,half,one,two
double precision b(m),sol(n2,nsol),a(m,nn),x(nn),wa(m5,n4),wb(m)
double precision sum,e(m),dsol(m,ndsol)
double precision qn(nn),cutoff,ci(4,nn),tnmat(4,nn),tnew,told,tn
data zero/0.00/
data half/0.50/
data one/1.00/
data two/2.00/
#
#  check dimension parameters
#
n=nn
ift = 0
wa(m+2,nn+1) = one
if (m5!=m+5)
	ift = 3
if (n2!=n+2)
	ift = 4
if (n4!=n+4)
	ift = 5
if (m<=zero||n<=zero)
	ift = 6
if (ift<=two) {
#
#  initialization
#
	iend = .true.
	lci2 = .false.
	lup = .true.
	skip = .false.
	idxcf = 0
	tnew = zero
	tn = zero
	m1 = m+1
	n1 = n+1
	n3 = n+3
	m2 = m+2
	m3 = m+3
	m4 = m+4
	do j = 1,n{
		x(j) = zero
		}
	do i = 1,m
		e(i) = zero
	if (t<zero||t>one) {
		t0 = one/float(m)-toler
		t1 = one-t0
		t = t0
		iend = .false.
		lci1 = .false.
		}
	repeat { #0
		do i = 1,m {
			k = 1
			do j = 1,nn
				if (k<=nn){
					if(j==idxcf)
						skip = .true.
					else
						skip = .false.
					if(!skip){
						wa(i,k) = a(i,j)
						k = k+1
						}
					}	
			wa(i,n4) = n+i
			wa(i,n2) = b(i)
			if (idxcf != 0)
				wa(i,n3) = tnew*a(i,idxcf)
			else
				wa(i,n3) = zero
			wa(i,n1) = wa(i,n2)-wa(i,n3)
			if (wa(i,n1)<zero)
				do j = 1,n4
					wa(i,j) = -wa(i,j)
			}
		do j = 1,n {
			wa(m4,j) = j
			wa(m2,j) = zero
			wa(m3,j) = zero
			do i = 1,m {
				aux = sign(one,wa(m4,j))*wa(i,j)
				wa(m2,j) = wa(m2,j)+aux*(one-sign(one,wa(i,n4)))
				wa(m3,j) = wa(m3,j)+aux*sign(one,wa(i,n4))
				}
			wa(m3,j) = two*wa(m3,j)
			}
		dif = zero
		init = .false.
		if(!lci2){
#compute the columns means
			do k = 1,n {
				wa(m5,k) = zero
				do i = 1,m{
					wa(m5,k) = wa(m5,k)+a(i,k)
					}
				wa(m5,k) = wa(m5,k)/float(m)
				}
			}
		lsol = 1
		kount = 0
		repeat { #1
#
# compute new marginal costs
#
			do j = 1,n{
				wa(m1,j) = wa(m2,j)+wa(m3,j)*t
				}
			if (!init) {
#
# stage 1
#
# determine the vector to enter the basis
#
				stage = .true.
				kr = 1
				kl = 1
				go to 30
				}
			repeat { #2
#
# stage 2
#
				stage = .false.
				repeat { #3
#
# determine the vector to enter the basis
#
					max = -big
					do j = kr,n {
						d = wa(m1,j)
						if (d<zero) {
							if (d>(-two))
								next 1
							d = -d-two
							}
						if (d>max) {
							max = d
							in = j
							}
						}
					if (max<=toler)
						break 2
					if (wa(m1,in)<=zero) {
						do i = 1,m4
							wa(i,in) = -wa(i,in)
						wa(m1,in) = wa(m1,in)-two
						wa(m2,in) = wa(m2,in)-two
						}
					repeat { #4
#
# determine the vector to leave the basis
#
						k = 0
						do i = kl,m {
							d = wa(i,in)
							if (d>toler) {
								k = k+1
								wb(k) = wa(i,n1)/d
								s(k) = i
								test = .true.
								}
							}
						repeat { #5
							if (k<=0)
								test = .false.
							else {
								min = big
								do i = 1,k
									if (wb(i)<min) {
										j = i
										min = wb(i)
										out = s(i)
										}
								wb(j) = wb(k)
								s(j) = s(k)
								k = k-1
								}
#
# check for linear dependence in stage 1
#
							if (!test&&stage)
								break 1
							if (!test)
								break 5
							pivot = wa(out,in)
							if (wa(m1,in)-pivot-pivot<=toler){
								go to 10
								}
							do j = kr,n3 {
								d = wa(out,j)
								wa(m1,j) = wa(m1,j)-d-d
								wa(m2,j) = wa(m2,j)-d-d
								wa(out,j) = -d
								}
							wa(out,n4) = -wa(out,n4)
							} #5
						do i = 1,m4 {
							d = wa(i,kr)
							wa(i,kr) = wa(i,in)
							wa(i,in) = d
							}
						kr = kr+1
						go to 20
#
# pivot on wa(out,in)
#
						10  do j = kr,n3
							if (j!=in)
								wa(out,j) = wa(out,j)/pivot
						do i = 1,m3
							if (i!=out) {
								d = wa(i,in)
								do j = kr,n3
									if (j!=in)
										wa(i,j) = wa(i,j)-d*wa(out,j)
								}
						do i = 1,m3
							if (i!=out)
								wa(i,in) = -wa(i,in)/pivot
						wa(out,in) = one/pivot
						d = wa(out,n4)
						wa(out,n4) = wa(m4,in)
						wa(m4,in) = d
						kount = kount+1
						if (!stage)
							break 1
#
# interchange rows in stage 1
#
						kl = kl+1
						do j = kr,n4 {
							d = wa(out,j)
							wa(out,j) = wa(kount,j)
							wa(kount,j) = d
							}
						20  if (kount+kr==n1)
							break 2
						30  max = -one
						do j = kr,n
							if (abs(wa(m4,j))<=n) {
								d = abs(wa(m1,j))
								if (d>max) {
									max = d
									in = j
									}
								}
						if (wa(m1,in)<zero)
							do i = 1,m4
								wa(i,in) = -wa(i,in)
						} #4
					} #3
				} #2
			if (kr==1) {
				do j = 1,n {
					d = abs(wa(m1,j))
					if (d<=toler||two-d<=toler){
						ift = 1
						wa(m2,nn+1) = zero
						go to 80
						}
					}
				}
			80 kount = 0
			sum = zero
			if (!lci2){
				do i = 1,kl-1 {
					k = wa(i,n4)*sign(one,wa(i,n4))
					x(k) = wa(i,n1)*sign(one,wa(i,n4))
					}
				}
			do i = 1,n {
				kd = abs(wa(m4,i))-n
				dsol(kd,lsol) = one+wa(m1,i)/two
				if (wa(m4,i)<zero)
					dsol(kd,lsol) = one-dsol(kd,lsol)
				if (!lci2){
					sum = sum+x(i)*wa(m5,i)
					sol(i+2,lsol) = x(i)
					h(i,lsol) = kd
					}
				}
			do i = kl,m {
				kd = abs(wa(i,n4))-n
				if (wa(i,n4)<zero)
					dsol(kd,lsol) = zero
				if (wa(i,n4)>zero)
					dsol(kd,lsol) = one
				}
			if (!lci2){
				sol(1,lsol) = smax
				sol(2,lsol) = sum
				do i=1,m
					dsol(i,lsol+1) = dsol(i,lsol)
				}
			if (lci2){
# compute next theta
				a1 = zero
				do i = 1,m {
					a1 = a1+a(i,idxcf)*(dsol(i,lsol)+t-one)
					}
				tn = a1/sqrt(qn(idxcf)*t*(one-t))
				if (abs(tn)<cutoff){
					if (lup)
						smax = big
					else
						smax = -big
					do i =1,kl-1 {
						k = wa(i,n4)*sign(one,wa(i,n4))
						sol(k,1) = wa(i,n2)*sign(one,wa(i,n4))
						sol(k,2) = wa(i,n3)*sign(one,wa(i,n4))/tnew
						}
					do i = kl,m {
						a1 = zero
						b1 = zero
						k = wa(i,n4)*sign(one,wa(i,n4))-n
						l = 1
						do j = 1,n{
							if (j==idxcf)
								l = l+1
							a1 = a1 + a(k,l)*sol(j,1)
							b1 = b1 + a(k,l)*sol(j,2)
							l = l+1
							}
						tnt = (b(k)-a1)/(a(k,idxcf)-b1)
						if (lup){
							if (tnt>tnew)
								if (tnt<smax){
									smax = tnt
									out = i
									}
							}
						else {
							if (tnt<tnew)
								if (tnt>smax){
									smax = tnt
									out = i
									}
							}
						}
					if (lup){
						told = tnew 
						tnew = smax + toler
						ci(3,idxcf) = told - toler
						tnmat(3,idxcf) = tn
						if (!(tnew < big-toler)){
							ci(3,idxcf) = big
							ci(4,idxcf) = big
							tnmat(3,idxcf) = tn
							tnmat(4,idxcf) = tn
							lup = .false.
							go to 70
							}
						}
					else{
						told = tnew 
						tnew = smax - toler
						ci(2,idxcf) = told + toler
						tnmat(2,idxcf) = tn
						if (!(tnew > -big+toler)){
							ci(2,idxcf) = -big
							ci(1,idxcf) = -big
							tnmat(2,idxcf) = tn
							tnmat(1,idxcf) = tn
							lup = .true.
							go to 60
							}
						}
#update the new marginal cost
					do i = 1,m{
						wa(i,n3) = wa(i,n3)/told*tnew
						wa(i,n1) = wa(i,n2) - wa(i,n3)
						}
					do j = kr,n3{
						d = wa(out,j)
						wa(m1,j) = wa(m1,j) -d -d
						wa(m2,j) = wa(m2,j) -d -d
						wa(out,j) = -d
						}
					wa(out,n4) = -wa(out,n4)
					init = .true.
					}
				else{
					if (lup){
						ci(4,idxcf) = tnew - toler
						tnmat(4,idxcf) = tn
						lup = .false.
						go to 70
						}
					else{
						ci(1,idxcf) = tnew + toler
						tnmat(1,idxcf) = tn
						lup = .true.
						go to 60
						}
					}
				}
			if ((iend)&&(!lci2))
				go to 40
			if (!lci2){
				init = .true.
				lsol = lsol+1
				do i = 1,m
					s(i) = zero
				do j = 1,n
					x(j) = zero
#
#  compute next t
#
				smax = two
				do j = 1,n {
					b1 = wa(m3,j)
					a1 = (-two-wa(m2,j))/b1
					b1 = -wa(m2,j)/b1
					if (a1>=t)
						if (a1<smax) {
							smax = a1
							dif = (b1-a1)/two
							}
					if (b1>t)
						if (b1<smax) {
							smax = b1
							dif = (b1-a1)/two
							}
					}
				tnt = smax+toler*(one+abs(dif))
				if (tnt>=t1+toler)
					iend = .true.
				t = tnt
				if (iend)
					t = t1
				}
			} #1
		wa(m2,nn+1) = two
		ift = 2
		go to 50
		40  if (lsol>2) {
			sol(1,1) = zero
			sol(1,lsol) = one
			do i = 1,m {
				dsol(i,1) = one
				dsol(i,lsol) = zero
				dsol(i,lsol+1) = zero
				}
			}
		l = kl-1
		do i = 1,l
			if (wa(i,n1)<zero)
				do j = kr,n4
					wa(i,j) = -wa(i,j)
		50  sum = zero
		if(!lci2){
			do i = kl,m {
				k = wa(i,n4)*sign(one,wa(i,n4))
				d = wa(i,n1)*sign(one,wa(i,n4))
				sum = sum+d*sign(one,d)*(half+sign(one,d)*(t-half))
				k = k-n
				e(k) = d
				}
			wa(m2,n2) = kount
			wa(m1,n2) = n1-kr
			wa(m1,n1) = sum
			}
		if (wa(m2,nn+1)==two)
			break 1
		if (!lci1)
			break 1
		if (!lci2){
			lci2 = .true.
			n = nn-1
			n1 = n+1
			n2 = n+2
			n3 = n+3
			n4 = n+4
		60	idxcf = idxcf+1
			if (idxcf>nn)
				break 1
		70	if (lup){
				tnew = x(idxcf)+toler
				told = tnew
				ci(3,idxcf) = x(idxcf)
				tnmat(3,idxcf) = zero
				}
			else{
				tnew = x(idxcf)-toler
				told = tnew
				ci(2,idxcf) = x(idxcf)
				tnmat(2,idxcf) = zero
				}
			}
		} #0
	# restore the original value of dual when ci is true
	do i=1,m
		dsol(i,lsol) = dsol(i,lsol+1)
	}
return
end