# This is a special version of K-dO rq modified to compute the LR process
# Note that the sol array is now p+3 by J and the third row contains the
# value of the objective function at the tau specified in row one.
#
#     This software is in the public domain  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 rqbr(m,nn,m5,n3,n4,a,b,t,toler,ift,x,e,s,wa,wb,nsol,ndsol,sol,dsol,lsol,h,qn,cutoff,ci,tnmat,big,lci1)
#
#     m = number of observations
#     n = number of parameters
#     m5 = m+5
#     n3 = n+3
#     n4 = 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 see below
#     x the parameter estimate betahat
#     e is the residual vector
#     s is an integer work array
#     wa is a real work array
#     wb is another real 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 confidence 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
#     neednt 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 dorey, 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,n3,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(n3,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
parameter( zero = 0.d0)
parameter( one = 1.d0)
parameter( two = 2.d0)
#
#  check dimension parameters
#
n=nn
ift = 0
wa(m+2,nn+1) = one
if (m5!=m+5)
  ift = 3
if (n3!=n+3)
  ift = 4
if (n4!=n+4)
  ift = 5
if (m<=zero||n<=zero)
  ift = 6
if (ift<=two) {
#
#  initialization
#
  half = one/two
  iend = .true.
  lci2 = .false.
  lup = .true.
  skip = .false.
  idxcf = 0
  tnew = zero
  tn = zero
  m1 = m+1
  n1 = n+1
  n2 = n+2
  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/(two*float(m))
    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+3,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
        sum = zero
        do j=kl,m{
          d = wa(j,n1)*sign(one,wa(j,n4))
          sum = sum + d*(smax + half*(sign(one,d) - one))
	 }
        sol(3,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(2,1) = zero
      sol(3,1) = zero
      sol(1,lsol) = one
      #sol(2,lsol) = zero
      sol(3,lsol) = zero
      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