      subroutine ssort(x,y,n,kflag)
c***begin prologue  ssort
c***date written   761101   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  n6a2b1
c***keywords  quicksort,singleton quicksort,sort,sorting
c***author  jones, r. e., (snla)
c           wisniewski, j. a., (snla)
c***purpose  ssort sorts array x and optionally makes the same
c            interchanges in array y.  the array x may be sorted in
c            increasing order or decreasing order.  a slightly modified
c            quicksort algorithm is used.
c***description
c
c     written by rondall e. jones
c     modified by john a. wisniewski to use the singleton quicksort
c     algorithm.  date 18 november 1976.
c
c     abstract
c         ssort sorts array x and optionally makes the same
c         interchanges in array y.  the array x may be sorted in
c         increasing order or decreasing order.  a slightly modified
c         quicksort algorithm is used.
c
c     reference
c         singleton, r. c., algorithm 347, an efficient algorithm for
c         sorting with minimal storage, cacm,12(3),1969,185-7.
c
c     description of parameters
c         x - array of values to be sorted   (usually abscissas)
c         y - array to be (optionally) carried along
c         n - number of values in array x to be sorted
c         kflag - control parameter
c             =2  means sort x in increasing order and carry y along.
c             =1  means sort x in increasing order (ignoring y)
c             =-1 means sort x in decreasing order (ignoring y)
c             =-2 means sort x in decreasing order and carry y along.
c***references  singleton,r.c., algorithm 347, an efficient algorithm
c                 for sorting with minimal storage, cacm,12(3),1969,
c                 185-7.
c***routines called  none
c***end prologue  ssort
      parameter (nerr=19)
      include 'parameters.h'
      dimension x(*),y(*),il(21),iu(21)
c***first executable statement  ssort
      nn = n
      if (nn.ge.1) go to 10
      if(inter)
     $  write(6,*)' the number of values to be sorted was not positive'
      write(nerr,*)
     $        ' the number of values to be sorted was not positive'
      return
   10 kk = iabs(kflag)
      if ((kk.eq.1).or.(kk.eq.2)) go to 15
      return
c
c alter array x to get decreasing order if needed
c
   15 if (kflag.ge.1) go to 30
      do 20 i=1,nn
   20 x(i) = -x(i)
   30 go to (100,200),kk
c
c sort x only
c
  100 continue
      m=1
      i=1
      j=nn
      r=.375
  110 if (i .eq. j) go to 155
      if (r .gt. .5898437) go to 120
      r=r+3.90625e-2
      go to 125
  120 r=r-.21875
  125 k=i
c                                  select a central element of the
c                                  array and save it in location t
      ij = i + ifix (float (j-i) * r)
      t=x(ij)
c                                  if first element of array is greater
c                                  than t, interchange with t
      if (x(i) .le. t) go to 130
      x(ij)=x(i)
      x(i)=t
      t=x(ij)
  130 l=j
c                                  if last element of array is less thans
c                                  t, interchange with t
      if (x(j) .ge. t) go to 140
      x(ij)=x(j)
      x(j)=t
      t=x(ij)
c                                  if first element of array is greater
c                                  than t, interchange with t
      if (x(i) .le. t) go to 140
      x(ij)=x(i)
      x(i)=t
      t=x(ij)
      go to 140
  135 tt=x(l)
      x(l)=x(k)
      x(k)=tt
c                                  find an element in the second half ofs
c                                  the array which is smaller than t
  140 l=l-1
      if (x(l) .gt. t) go to 140
c                                  find an element in the first half of
c                                  the array which is greater than t
  145 k=k+1
      if (x(k) .lt. t) go to 145
c                                  interchange these elements
      if (k .le. l) go to 135
c                                  save upper and lower subscripts of
c                                  the array yet to be sorted
      if (l-i .le. j-k) go to 150
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 160
  150 il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 160
c                                  begin again on another portion of
c                                  the unsorted array
  155 m=m-1
      if (m .eq. 0) go to 300
      i=il(m)
      j=iu(m)
  160 if (j-i .ge. 1) go to 125
      if (i .eq. 1) go to 110
      i=i-1
  165 i=i+1
      if (i .eq. j) go to 155
      t=x(i+1)
      if (x(i) .le. t) go to 165
      k=i
  170 x(k+1)=x(k)
      k=k-1
      if (t .lt. x(k)) go to 170
      x(k+1)=t
      go to 165
c
c sort x and carry y along
c
  200 continue
      m=1
      i=1
      j=nn
      r=.375
  210 if (i .eq. j) go to 255
      if (r .gt. .5898437) go to 220
      r=r+3.90625e-2
      go to 225
  220 r=r-.21875
  225 k=i
c                                  select a central element of the
c                                  array and save it in location t
      ij = i + ifix (float (j-i) *r)
      t=x(ij)
      ty= y(ij)
c                                  if first element of array is greater
c                                  than t, interchange with t
      if (x(i) .le. t) go to 230
      x(ij)=x(i)
      x(i)=t
      t=x(ij)
       y(ij)= y(i)
       y(i)=ty
      ty= y(ij)
  230 l=j
c                                  if last element of array is less thans
c                                  t, interchange with t
      if (x(j) .ge. t) go to 240
      x(ij)=x(j)
      x(j)=t
      t=x(ij)
       y(ij)= y(j)
       y(j)=ty
      ty= y(ij)
c                                  if first element of array is greater
c                                  than t, interchange with t
      if (x(i) .le. t) go to 240
      x(ij)=x(i)
      x(i)=t
      t=x(ij)
       y(ij)= y(i)
       y(i)=ty
      ty= y(ij)
      go to 240
  235 tt=x(l)
      x(l)=x(k)
      x(k)=tt
      tty= y(l)
       y(l)= y(k)
       y(k)=tty
c                                  find an element in the second half ofs
c                                  the array which is smaller than t
  240 l=l-1
      if (x(l) .gt. t) go to 240
c                                  find an element in the first half of
c                                  the array which is greater than t
  245 k=k+1
      if (x(k) .lt. t) go to 245
c                                  interchange these elements
      if (k .le. l) go to 235
c                                  save upper and lower subscripts of
c                                  the array yet to be sorted
      if (l-i .le. j-k) go to 250
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 260
  250 il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 260
c                                  begin again on another portion of
c                                  the unsorted array
  255 m=m-1
      if (m .eq. 0) go to 300
      i=il(m)
      j=iu(m)
  260 if (j-i .ge. 1) go to 225
      if (i .eq. 1) go to 210
      i=i-1
  265 i=i+1
      if (i .eq. j) go to 255
      t=x(i+1)
      ty= y(i+1)
      if (x(i) .le. t) go to 265
      k=i
  270 x(k+1)=x(k)
       y(k+1)= y(k)
      k=k-1
      if (t .lt. x(k)) go to 270
      x(k+1)=t
       y(k+1)=ty
      go to 265
c
c clean up
c
  300 if (kflag.ge.1) return
      do 310 i=1,nn
  310 x(i) = -x(i)
      return
      end
      subroutine dsort (dx, dy, n, kflag)
c***begin prologue  dsort
c***purpose  sort an array and optionally make the same interchanges in
c            an auxiliary array.  the array may be sorted in increasing
c            or decreasing order.  a slightly modified quicksort
c            algorithm is used.
c***library   slatec
c***category  n6a2b
c***type      double precision (ssort-s, dsort-d, isort-i)
c***keywords  singleton quicksort, sort, sorting
c***author  jones, r. e., (snla)
c           wisniewski, j. a., (snla)
c***description
c
c   dsort sorts array dx and optionally makes the same interchanges in
c   array dy.  the array dx may be sorted in increasing order or
c   decreasing order.  a slightly modified quicksort algorithm is used.
c
c   description of parameters
c      dx - array of values to be sorted   (usually abscissas)
c      dy - array to be (optionally) carried along
c      n  - number of values in array dx to be sorted
c      kflag - control parameter
c            =  2  means sort dx in increasing order and carry dy along.
c            =  1  means sort dx in increasing order (ignoring dy)
c            = -1  means sort dx in decreasing order (ignoring dy)
c            = -2  means sort dx in decreasing order and carry dy along.
c
c***references  r. c. singleton, algorithm 347, an efficient algorithm
c                 for sorting with minimal storage, communications of
c                 the acm, 12, 3 (1969), pp. 185-187.
c***end prologue  dsort
      parameter (nerr=19)
      include 'parameters.h'
      double precision dx(*), dy(*)
      double precision r, t, tt, tty, ty
      integer il(21), iu(21)
c***first executable statement  dsort
      nn = n
      if(nn.ge.1)goto 5
      if(inter)
     $  write(6,*)' the number of values to be sorted was not positive'
      write(nerr,*)
     $        ' the number of values to be sorted was not positive'
      return
c
    5 kk = abs(kflag)
      if((kk.eq.1).or.(kk.eq.2))goto 7
      return
c
c     alter array dx to get decreasing order if needed
c
    7 if (kflag .le. -1) then
         do 10 i=1,nn
            dx(i) = -dx(i)
   10    continue
      endif
c
      if (kk .eq. 2) go to 100
c
c     sort dx only
c
      m = 1
      i = 1
      j = nn
      r = 0.375d0
c
   20 if (i .eq. j) go to 60
      if (r .le. 0.5898437d0) then
         r = r+3.90625d-2
      else
         r = r-0.21875d0
      endif
c
   30 k = i
c
c     select a central element of the array and save it in location t
c
      ij = i + int((j-i)*r)
      t = dx(ij)
c
c     if first element of array is greater than t, interchange with t
c
      if (dx(i) .gt. t) then
         dx(ij) = dx(i)
         dx(i) = t
         t = dx(ij)
      endif
      l = j
c
c     if last element of array is less than than t, interchange with t
c
      if (dx(j) .lt. t) then
         dx(ij) = dx(j)
         dx(j) = t
         t = dx(ij)
c
c        if first element of array is greater than t, interchange with t
c
         if (dx(i) .gt. t) then
            dx(ij) = dx(i)
            dx(i) = t
            t = dx(ij)
         endif
      endif
c
c     find an element in the second half of the array which is smaller
c     than t
c
   40 l = l-1
      if (dx(l) .gt. t) go to 40
c
c     find an element in the first half of the array which is greater
c     than t
c
   50 k = k+1
      if (dx(k) .lt. t) go to 50
c
c     interchange these elements
c
      if (k .le. l) then
         tt = dx(l)
         dx(l) = dx(k)
         dx(k) = tt
         go to 40
      endif
c
c     save upper and lower subscripts of the array yet to be sorted
c
      if (l-i .gt. j-k) then
         il(m) = i
         iu(m) = l
         i = k
         m = m+1
      else
         il(m) = k
         iu(m) = j
         j = l
         m = m+1
      endif
      go to 70
c
c     begin again on another portion of the unsorted array
c
   60 m = m-1
      if (m .eq. 0) go to 190
      i = il(m)
      j = iu(m)
c
   70 if (j-i .ge. 1) go to 30
      if (i .eq. 1) go to 20
      i = i-1
c
   80 i = i+1
      if (i .eq. j) go to 60
      t = dx(i+1)
      if (dx(i) .le. t) go to 80
      k = i
c
   90 dx(k+1) = dx(k)
      k = k-1
      if (t .lt. dx(k)) go to 90
      dx(k+1) = t
      go to 80
c
c     sort dx and carry dy along
c
  100 m = 1
      i = 1
      j = nn
      r = 0.375d0
c
  110 if (i .eq. j) go to 150
      if (r .le. 0.5898437d0) then
         r = r+3.90625d-2
      else
         r = r-0.21875d0
      endif
c
  120 k = i
c
c     select a central element of the array and save it in location t
c
      ij = i + int((j-i)*r)
      t = dx(ij)
      ty = dy(ij)
c
c     if first element of array is greater than t, interchange with t
c
      if (dx(i) .gt. t) then
         dx(ij) = dx(i)
         dx(i) = t
         t = dx(ij)
         dy(ij) = dy(i)
         dy(i) = ty
         ty = dy(ij)
      endif
      l = j
c
c     if last element of array is less than t, interchange with t
c
      if (dx(j) .lt. t) then
         dx(ij) = dx(j)
         dx(j) = t
         t = dx(ij)
         dy(ij) = dy(j)
         dy(j) = ty
         ty = dy(ij)
c
c        if first element of array is greater than t, interchange with t
c
         if (dx(i) .gt. t) then
            dx(ij) = dx(i)
            dx(i) = t
            t = dx(ij)
            dy(ij) = dy(i)
            dy(i) = ty
            ty = dy(ij)
         endif
      endif
c
c     find an element in the second half of the array which is smaller
c     than t
c
  130 l = l-1
      if (dx(l) .gt. t) go to 130
c
c     find an element in the first half of the array which is greater
c     than t
c
  140 k = k+1
      if (dx(k) .lt. t) go to 140
c
c     interchange these elements
c
      if (k .le. l) then
         tt = dx(l)
         dx(l) = dx(k)
         dx(k) = tt
         tty = dy(l)
         dy(l) = dy(k)
         dy(k) = tty
         go to 130
      endif
c
c     save upper and lower subscripts of the array yet to be sorted
c
      if (l-i .gt. j-k) then
         il(m) = i
         iu(m) = l
         i = k
         m = m+1
      else
         il(m) = k
         iu(m) = j
         j = l
         m = m+1
      endif
      go to 160
c
c     begin again on another portion of the unsorted array
c
  150 m = m-1
      if (m .eq. 0) go to 190
      i = il(m)
      j = iu(m)
c
  160 if (j-i .ge. 1) go to 120
      if (i .eq. 1) go to 110
      i = i-1
c
  170 i = i+1
      if (i .eq. j) go to 150
      t = dx(i+1)
      ty = dy(i+1)
      if (dx(i) .le. t) go to 170
      k = i
c
  180 dx(k+1) = dx(k)
      dy(k+1) = dy(k)
      k = k-1
      if (t .lt. dx(k)) go to 180
      dx(k+1) = t
      dy(k+1) = ty
      go to 170
c
c     clean up
c
  190 if (kflag .le. -1) then
         do 200 i=1,nn
            dx(i) = -dx(i)
  200    continue
      endif
      return
      end
      subroutine fzero(f,b,c,r,re,ae,iflag,rpar,ipar)
c***begin prologue  fzero
c***date written   700901   (yymmdd)
c***revision date  820801   (yymmdd)
c***category no.  f1b
c***keywords  bisection,nonlinear,roots,zeros
c***author  shampine,l.f.,snla
c           watts,h.a.,snla
c***purpose  fzero searches for a zero of a function f(x) in a given
c            interval (b,c).  it is designed primarily for problems
c            where f(b) and f(c) have opposite signs.
c***description
c
c     based on a method by t j dekker
c     written by l f shampine and h a watts
c
c            fzero searches for a zero of a function f(x) between
c            the given values b and c until the width of the interval
c            (b,c) has collapsed to within a tolerance specified by
c            the stopping criterion, abs(b-c) .le. 2.*(rw*abs(b)+ae).
c            the method used is an efficient combination of bisection
c            and the secant rule.
c
c     description of arguments
c
c     f,b,c,r,re and ae are input parameters
c     b,c and iflag are output parameters (flagged by an * below)
c
c        f     - name of the real valued external function.  this name
c                must be in an external statement in the calling
c                program.  f must be a function of one real argument.
c
c       *b     - one end of the interval (b,c).  the value returned for
c                b usually is the better approximation to a zero of f.
c
c       *c     - the other end of the interval (b,c)
c
c        r     - a (better) guess of a zero of f which could help in
c                speeding up convergence.  if f(b) and f(r) have
c                opposite signs, a root will be found in the interval
c                (b,r); if not, but f(r) and f(c) have opposite
c                signs, a root will be found in the interval (r,c);
c                otherwise, the interval (b,c) will be searched for a
c                possible root.  when no better guess is known, it is
c                recommended that r be set to b or c; because if r is
c                not interior to the interval (b,c), it will be ignored.
c
c        re    - relative error used for rw in the stopping criterion.
c                if the requested re is less than machine precision,
c                then rw is set to approximately machine precision.
c
c        ae    - absolute error used in the stopping criterion.  if the
c                given interval (b,c) contains the origin, then a
c                nonzero value should be chosen for ae.
c
c       *iflag - a status code.  user must check iflag after each call.
c                control returns to the user from fzero in all cases.
c                xerror does not process diagnostics in these cases.
c
c                1  b is within the requested tolerance of a zero.
c                   the interval (b,c) collapsed to the requested
c                   tolerance, the function changes sign in (b,c), and
c                   f(x) decreased in magnitude as (b,c) collapsed.
c
c                2  f(b) = 0.  however, the interval (b,c) may not have
c                   collapsed to the requested tolerance.
c
c                3  b may be near a singular point of f(x).
c                   the interval (b,c) collapsed to the requested tol-
c                   erance and the function changes sign in (b,c), but
c                   f(x) increased in magnitude as (b,c) collapsed,i.e.
c                     abs(f(b out)) .gt. max(abs(f(b in)),abs(f(c in)))
c
c                4  no change in sign of f(x) was found although the
c                   interval (b,c) collapsed to the requested tolerance.
c                   the user must examine this case and decide whether
c                   b is near a local minimum of f(x), or b is near a
c                   zero of even multiplicity, or neither of these.
c
c                5  too many (.gt. 500) function evaluations used.
c      *ipar     vector of integer parameters to pass to f
c***references  l. f. shampine and h. a. watts, *fzero, a root-solving
c                 code*, sc-tm-70-631, september 1970.
c               t. j. dekker, *finding a zero by means of successive
c                 linear interpolation*, 'constructive aspects of the
c                 fundamental theorem of algebra', edited by b. dejon
c                 p. henrici, 1969.
c***routines called  r1mach
c***end prologue  fzero
c
      dimension rpar(*)
      dimension ipar(*)
c
c     er is two times the computer unit roundoff value which is
c     defined here by the function r1mach.

c***first executable statement  fzero
      er = 2.0 * r1mach(4)
c
c     initialize
c
      z=r
      if(r.le.min(b,c).or.r.ge.max(b,c)) z=c
      rw=max(re,er)
      aw=max(ae,0.)
      ic=0
      t=z
      fz=f(t,rpar,ipar)
      t=b
      fb=f(t,rpar,ipar)
      kount=2
      if(sign(1.0,fz).eq.sign(1.0,fb)) go to 1
      c=z
      fc=fz
      go to 2
    1 if(z.eq.c) go to 2
      t=c
      fc=f(t,rpar,ipar)
      kount=3
      if(sign(1.0,fz).eq.sign(1.0,fc)) go to 2
      b=z
      fb=fz
    2 a=c
      fa=fc
      acbs=abs(b-c)
      fx=max(abs(fb),abs(fc))
c
    3 if (abs(fc) .ge. abs(fb)) go to 4
c     perform interchange
      a=b
      fa=fb
      b=c
      fb=fc
      c=a
      fc=fa
c
    4 cmb=0.5*(c-b)
      acmb=abs(cmb)
      tol=rw*abs(b)+aw
c
c     test stopping criterion and function count
c
      if (acmb .le. tol) go to 10
      if(fb.eq.0.) go to 11
      if(kount.ge.500) go to 14
c
c     calculate new iterate implicitly as b+p/q
c     where we arrange p .ge. 0.
c     the implicit form is used to prevent overflow.
c
      p=(b-a)*fb
      q=fa-fb
      if (p .ge. 0.) go to 5
      p=-p
      q=-q
c
c     update a and check for satisfactory reduction
c     in the size of the bracketing interval.
c     if not, perform bisection.
c
    5 a=b
      fa=fb
      ic=ic+1
      if (ic .lt. 4) go to 6
      if (8.*acmb .ge. acbs) go to 8
      ic=0
      acbs=acmb
c
c     test for too small a change
c
    6 if (p .gt. abs(q)*tol) go to 7
c
c     increment by tolerance
c
      b=b+sign(tol,cmb)
      go to 9
c
c     root ought to be between b and (c+b)/2.
c
    7 if (p .ge. cmb*q) go to 8
c
c     use secant rule
c
      b=b+p/q
      go to 9
c
c     use bisection
c
    8 b=0.5*(c+b)
c
c     have completed computation for new iterate b
c
    9 t=b
      fb=f(t,rpar,ipar)
      kount=kount+1
c
c     decide whether next step is interpolation or extrapolation
c
      if (sign(1.0,fb) .ne. sign(1.0,fc)) go to 3
      c=a
      fc=fa
      go to 3
c
c
c     finished. process results for proper setting of iflag
c
   10 if (sign(1.0,fb) .eq. sign(1.0,fc)) go to 13
      if (abs(fb) .gt. fx) go to 12
      iflag = 1
      return
   11 iflag = 2
      return
   12 iflag = 3
      return
   13 iflag = 4
      return
   14 iflag = 5
      return
      end
      subroutine dpqr79 (ndeg, coeff, root, ierr, work)
c***purpose  find the zeros of a polynomial with real coefficients.
c***library   slatec
c***type      single precision (rpqr79-s, cpqr79-c), converted to double precisi
c***description
c
c   abstract
c       this routine computes all zeros of a polynomial of degree ndeg
c       with real coefficients by computing the eigenvalues of the
c       companion matrix.
c
c   description of parameters
c       the user must dimension all arrays appearing in the call list
c            coeff(ndeg+1), root(ndeg), work(ndeg*(ndeg+2))
c
c    --input--
c      ndeg    degree of polynomial
c
c      coeff   real coefficients in descending order.  i.e.,
c              p(z)= coeff(1)*(z**ndeg) + coeff(ndeg)*z + coeff(ndeg+1)
c
c      work    real work array of dimension at least ndeg*(ndeg+2)
c
c   --output--
c      root    complex vector of roots
c
c      ierr    output error code
c           - normal code
c          0  means the roots were computed.
c           - abnormal codes
c          1  more than 30 qr iterations on some eigenvalue of the
c             companion matrix
c          2  coeff(1)=0.0
c          3  ndeg is invalid (less than or equal to 0)
c
      double precision coeff(*), work(*), scale
      double complex root(*)
      integer ndeg, ierr, k, kh, kwr, kwi, kcol
      ierr = 0
      if (abs(coeff(1)) .eq. 0.0) then
         ierr = 2
         return
      endif
c
      if (ndeg .le. 0) then
         ierr = 3
         return
      endif
c
      if (ndeg .eq. 1) then
         root(1) = cmplx(-coeff(2)/coeff(1),0.0)
         return
      endif
c
      scale = 1.0e0/coeff(1)
      kh = 1
      kwr = kh+ndeg*ndeg
      kwi = kwr+ndeg
      kwend = kwi+ndeg-1
c
      do 10 k=1,kwend
         work(k) = 0.0e0
   10 continue
c
      do 20 k=1,ndeg
         kcol = (k-1)*ndeg+1
         work(kcol) = -coeff(k+1)*scale
         if (k .ne. ndeg) work(kcol+k) = 1.0e0
   20 continue
c
      call hqr (ndeg,ndeg,1,ndeg,work(kh),work(kwr),work(kwi),ierr)
c
      if (ierr .ne. 0) then
         ierr = 1
         return
      endif
c
      do 30 k=1,ndeg
         km1 = k-1
         root(k) = cmplx(work(kwr+km1),work(kwi+km1))
   30 continue
      return
      end
      subroutine hqr (nm, n, low, igh, h, wr, wi, ierr)
c***purpose  compute the eigenvalues of a real upper hessenberg matrix
c            using the qr method.
c***library   slatec (eispack)
c
c     this subroutine is a translation of the algol procedure hqr,
c     num. math. 14, 219-231(1970) by martin, peters, and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 359-371(1971).
c
c     this subroutine finds the eigenvalues of a real
c     upper hessenberg matrix by the qr method.
c
c     on input
c
c        nm must be set to the row dimension of the two-dimensional
c          array parameter, h, as declared in the calling program
c          dimension statement.  nm is an integer variable.
c
c        n is the order of the matrix h.  n is an integer variable.
c          n must be less than or equal to nm.
c
c        low and igh are two integer variables determined by the
c          balancing subroutine  balanc.  if  balanc  has not been
c          used, set low=1 and igh equal to the order of the matrix, n.
c
c        h contains the upper hessenberg matrix.  information about
c          the transformations used in the reduction to hessenberg
c          form by  elmhes  or  orthes, if performed, is stored
c          in the remaining triangle under the hessenberg matrix.
c          h is a two-dimensional real array, dimensioned h(nm,n).
c
c     on output
c
c        h has been destroyed.  therefore, it must be saved before
c          calling  hqr  if subsequent calculation and back
c          transformation of eigenvectors is to be performed.
c
c        wr and wi contain the real and imaginary parts, respectively,
c          of the eigenvalues.  the eigenvalues are unordered except
c          that complex conjugate pairs of values appear consecutively
c          with the eigenvalue having the positive imaginary part first.
c          if an error exit is made, the eigenvalues should be correct
c          for indices ierr+1, ierr+2, ..., n.  wr and wi are one-
c          dimensional real arrays, dimensioned wr(n) and wi(n).
c
c        ierr is an integer flag set to
c          zero       for normal return,
c          j          if the j-th eigenvalue has not been
c                     determined after a total of 30*n iterations.
c                     the eigenvalues should be correct for indices
c                     ierr+1, ierr+2, ..., n.
c
c     ------------------------------------------------------------------
c
c
      integer i,j,k,l,m,n,en,ll,mm,na,nm,igh,itn,its,low,mp2,enm2,ierr
      double precision h(nm,*),wr(*),wi(*)
      double precision p,q,r,s,t,w,x,y,zz,norm,s1,s2
      logical notlas
c
c***first executable statement  hqr
      ierr = 0
      norm = 0.0e0
      k = 1
c     .......... store roots isolated by balanc
c                and compute matrix norm ..........
      do 50 i = 1, n
c
         do 40 j = k, n
   40    norm = norm + abs(h(i,j))
c
         k = i
         if (i .ge. low .and. i .le. igh) go to 50
         wr(i) = h(i,i)
         wi(i) = 0.0e0
   50 continue
c
      en = igh
      t = 0.0e0
      itn = 30*n
c     .......... search for next eigenvalues ..........
   60 if (en .lt. low) go to 1001
      its = 0
      na = en - 1
      enm2 = na - 1
c     .......... look for single small sub-diagonal element
c                for l=en step -1 until low do -- ..........
   70 do 80 ll = low, en
         l = en + low - ll
         if (l .eq. low) go to 100
         s = abs(h(l-1,l-1)) + abs(h(l,l))
         if (s .eq. 0.0e0) s = norm
         s2 = s + abs(h(l,l-1))
         if (s2 .eq. s) go to 100
   80 continue
c     .......... form shift ..........
  100 x = h(en,en)
      if (l .eq. en) go to 270
      y = h(na,na)
      w = h(en,na) * h(na,en)
      if (l .eq. na) go to 280
      if (itn .eq. 0) go to 1000
      if (its .ne. 10 .and. its .ne. 20) go to 130
c     .......... form exceptional shift ..........
      t = t + x
c
      do 120 i = low, en
  120 h(i,i) = h(i,i) - x
c
      s = abs(h(en,na)) + abs(h(na,enm2))
      x = 0.75e0 * s
      y = x
      w = -0.4375e0 * s * s
  130 its = its + 1
      itn = itn - 1
c     .......... look for two consecutive small
c                sub-diagonal elements.
c                for m=en-2 step -1 until l do -- ..........
      do 140 mm = l, enm2
         m = enm2 + l - mm
         zz = h(m,m)
         r = x - zz
         s = y - zz
         p = (r * s - w) / h(m+1,m) + h(m,m+1)
         q = h(m+1,m+1) - zz - r - s
         r = h(m+2,m+1)
         s = abs(p) + abs(q) + abs(r)
         p = p / s
         q = q / s
         r = r / s
         if (m .eq. l) go to 150
         s1 = abs(p) * (abs(h(m-1,m-1)) + abs(zz) + abs(h(m+1,m+1)))
         s2 = s1 + abs(h(m,m-1)) * (abs(q) + abs(r))
         if (s2 .eq. s1) go to 150
  140 continue
c
  150 mp2 = m + 2
c
      do 160 i = mp2, en
         h(i,i-2) = 0.0e0
         if (i .eq. mp2) go to 160
         h(i,i-3) = 0.0e0
  160 continue
c     .......... double qr step involving rows l to en and
c                columns m to en ..........
      do 260 k = m, na
         notlas = k .ne. na
         if (k .eq. m) go to 170
         p = h(k,k-1)
         q = h(k+1,k-1)
         r = 0.0e0
         if (notlas) r = h(k+2,k-1)
         x = abs(p) + abs(q) + abs(r)
         if (x .eq. 0.0e0) go to 260
         p = p / x
         q = q / x
         r = r / x
  170    s = sign(sqrt(p*p+q*q+r*r),p)
         if (k .eq. m) go to 180
         h(k,k-1) = -s * x
         go to 190
  180    if (l .ne. m) h(k,k-1) = -h(k,k-1)
  190    p = p + s
         x = p / s
         y = q / s
         zz = r / s
         q = q / p
         r = r / p
c     .......... row modification ..........
         do 210 j = k, en
            p = h(k,j) + q * h(k+1,j)
            if (.not. notlas) go to 200
            p = p + r * h(k+2,j)
            h(k+2,j) = h(k+2,j) - p * zz
  200       h(k+1,j) = h(k+1,j) - p * y
            h(k,j) = h(k,j) - p * x
  210    continue
c
         j = min(en,k+3)
c     .......... column modification ..........
         do 230 i = l, j
            p = x * h(i,k) + y * h(i,k+1)
            if (.not. notlas) go to 220
            p = p + zz * h(i,k+2)
            h(i,k+2) = h(i,k+2) - p * r
  220       h(i,k+1) = h(i,k+1) - p * q
            h(i,k) = h(i,k) - p
  230    continue
c
  260 continue
c
      go to 70
c     .......... one root found ..........
  270 wr(en) = x + t
      wi(en) = 0.0e0
      en = na
      go to 60
c     .......... two roots found ..........
  280 p = (y - x) / 2.0e0
      q = p * p + w
      zz = sqrt(abs(q))
      x = x + t
      if (q .lt. 0.0e0) go to 320
c     .......... real pair ..........
      zz = p + sign(zz,p)
      wr(na) = x + zz
      wr(en) = wr(na)
      if (zz .ne. 0.0e0) wr(en) = x - w / zz
      wi(na) = 0.0e0
      wi(en) = 0.0e0
      go to 330
c     .......... complex pair ..........
  320 wr(na) = x + p
      wr(en) = x + p
      wi(na) = zz
      wi(en) = -zz
  330 en = enm2
      go to 60
c     .......... set error -- no convergence to an
c                eigenvalue after 30*n iterations ..........
 1000 ierr = en
 1001 return
      end
      function qbeta(pr,p,q)
      integer p,q,ipar(2)
      dimension rpar(1)
      external fbeta
      data re/1.e-6/,ae/1.e-30/
      xq=q
      if(p.eq.1)then
        qbeta=1.-(1.-pr)**(1./xq)
        return
      else
        xp=p
        ipar(1)=p
        ipar(2)=q
        rpar(1)=pr
        b=0.
        r=2.*xp/(xp+xq)
        c=5.*r
        call fzero(fbeta,b,c,r,re,ae,iflag,rpar,ipar)
        qbeta=b
        return
      endif
      end
      function fbeta(x,rpar,ipar)
      integer p,q
      double precision xx,xp,xq,xnum,xden,scale,sum,term
      dimension rpar(1),ipar(2)
      p=ipar(1)
      q=ipar(2)
      xx=x
      xp=p
      xq=q
      xnum=xq+xp-1.d0
      do 10 i=2,p 
   10   xnum=xnum*(xq+xp-i)
      xden=xp-1.d0
      do 11 i=p-2,2,-1
   11   xden=xden*i
      scale=xnum/xden
      sum=0.d0
      do 30 i=1,p
        xnum=1.d0
        do 20 j=1,i-1
   20     xnum=xnum*(xp-j)
        xden=xq
        do 25 j=1,i-1
   25     xden=xden*(xq+j) 
        term=xnum/xden*xx**(p-i)*(1.d0-xx)**(q+i-1)
   30   sum=sum+term
      fbeta=1.d0-scale*sum-rpar(1)
      return
      end
      function betai (x, pin, qin)
c august 1980 version.  w. fullerton, c3, los alamos scientific lab.
c based on bosten and battiste, remark on algorithm 179, comm. acm,
c v 17, p 153, (1974).
c
c             input arguments --
c x      upper limit of integration.  x must be in (0,1) inclusive.
c p      first beta distribution parameter.  p must be gt 0.0.
c q      second beta distribution parameter.  q must be gt 0.0.
c betai  the incomplete beta function ratio is the probability that a
c        random variable from a beta distribution having parameters
c        p and q will be less than or equal to x.
c
      external albeta, alog, r1mach
      data eps, alneps, sml, alnsml / 4*0.0 /
c
      if (eps.ne.0.) go to 10
      eps = r1mach(3)
      alneps = alog(eps)
      sml = r1mach(1)
      alnsml = alog(sml)
c
 10   continue
c      if (x.lt.0. .or. x.gt.1.0) call seteru (
c     1  35hbetai   x is not in the range (0,1), 35, 1, 2)
c      if (pin.le.0. .or. qin.le.0.) call seteru (
c     1  29hbetai   a and/or b is le zero, 29, 2, 2)
c
      y = x
      p = pin
      q = qin
      if (q.le.p .and. x.lt.0.8) go to 20
      if (x.lt.0.2) go to 20
      y = 1.0 - y
      p = qin
      q = pin
c
 20   if ((p+q)*y/(p+1.).lt.eps) go to 80
c
c evaluate the infinite sum first.
c term will equal y**p/beta(ps,p) * (1.-ps)i * y**i / fac(i)
c
      ps = q - aint(q)
      if (ps.eq.0.) ps = 1.0
      xb = p*alog(y) -  albeta(ps, p) - alog(p)
      betai = 0.0
      if (xb.lt.alnsml) go to 40
c
      betai = exp (xb)
      term = betai*p
      if (ps.eq.1.0) go to 40
c
      n = amax1 (alneps/alog(y), 4.0)
      do 30 i=1,n
        term = term*(float(i)-ps)*y/float(i)
        betai = betai + term/(p+float(i))
 30   continue
c
c now evaluate the finite sum, maybe.
c
 40   if (q.le.1.0) go to 70
c
      xb = p*alog(y) + q*alog(1.0-y) - albeta(p,q) - alog(q)
      ib = amax1 (xb/alnsml, 0.0)
      term = exp (xb - float(ib)*alnsml)
      c = 1.0/(1.0-y)
      p1 = q*c/(p+q-1.)
c
      finsum = 0.0
      n = q
      if (q.eq.float(n)) n = n - 1
      do 50 i=1,n
        if (p1.le.1.0 .and. term/eps.le.finsum) go to 60
        term = (q-float(i-1))*c*term/(p+q-float(i))
c
        if (term.gt.1.0) ib = ib - 1
        if (term.gt.1.0) term = term*sml
c
        if (ib.eq.0) finsum = finsum + term
 50   continue
c
 60   betai = betai + finsum
 70   if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
      betai = amax1 (amin1 (betai, 1.0), 0.0)
      return
c
 80   betai = 0.0
      xb = p*alog(amax1(y,sml)) - alog(p) - albeta(p,q)
      if (xb.gt.alnsml .and. y.ne.0.) betai = exp (xb)
      if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
      return
c
      end
      function albeta (a, b)
c july 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      external alngam, alnrel, alog, gamma, r9lgmc
      data sq2pil / 0.9189385332 0467274 e0 /
c
      p = amin1 (a, b)
      q = amax1 (a, b)
c
c      if (p.le.0.0) call seteru (
c     1  38halbeta  both arguments must be gt zero, 38, 1, 2)
      if (p.ge.10.0) go to 30
      if (q.ge.10.0) go to 20
c
c p and q are small.
c
      albeta = alog(gamma(p) * (gamma(q)/gamma(p+q)) )
      return
c
c p is small, but q is big.
c
 20   corr = r9lgmc(q) - r9lgmc(p+q)
      albeta = alngam(p) + corr + p - p*alog(p+q) +
     1  (q-0.5)*alnrel(-p/(p+q))
      return
c
c p and q are big.
c
 30   corr = r9lgmc(p) + r9lgmc(q) - r9lgmc(p+q)
      albeta = -0.5*alog(q) + sq2pil + corr + (p-0.5)*alog(p/(p+q))
     1  + q*alnrel(-p/(p+q))
      return
c
      end
      function alog (x)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension alncs(6), center(4), alncen(5)
      external csevl, inits, r1mach
c
c series for aln        on the interval  0.          to  3.46021d-03
c                                        with weighted error   1.50e-16
c                                         log weighted error  15.82
c                               significant figures required  15.65
c                                    decimal places required  16.21
c
      data aln cs( 1) /   1.3347199877 973882e0 /
      data aln cs( 2) /    .0006937562 83284112e0 /
      data aln cs( 3) /    .0000004293 40390204e0 /
      data aln cs( 4) /    .0000000002 89338477e0 /
      data aln cs( 5) /    .0000000000 00205125e0 /
      data aln cs( 6) /    .0000000000 00000150e0 /
c
      data center(1) / 1.0 /
      data center(2) / 1.25 /
      data center(3) / 1.50 /
      data center(4) / 1.75 /
c
      data alncen(  1) / 0.0e0                                         /
      data alncen(  2) / +.2231435513 14209755 e+0                     /
      data alncen(  3) / +.4054651081 08164381 e+0                     /
      data alncen(  4) / +.5596157879 35422686 e+0                     /
      data alncen(  5) / +.6931471805 59945309 e+0                     /
c
c aln2 = alog(2.0) - 0.625
      data aln2 / 0.0681471805 59945309e0 /
      data nterms / 0 /
c
      if (nterms.eq.0) nterms = inits (alncs, 6, 28.9*r1mach(3))
c
c      if (x.le.0.) call seteru (
c     1  29halog    x is zero or negative, 29, 1, 2)
c
      call r9upak (x, y, n)
c
      xn = n - 1
      y = 2.0*y
      ntrval = 4.0*y - 2.5
      if (ntrval.eq.5) t = ((y-1.0)-1.0) / (y+2.0)
      if (ntrval.lt.5) t = (y-center(ntrval))/(y+center(ntrval))
      t2 = t*t
c
      alog = 0.625*xn + (aln2*xn + alncen(ntrval) + 2.0*t +
     1  t*t2*csevl(578.0*t2-1.0, alncs, nterms) )
c
      return
      end
      function alngam (x)
c august 1980 edition.   w. fullerton, c3, los alamos scientific lab.
c
      external alog, gamma, r1mach, r9lgmc
      data sq2pil / 0.9189385332 0467274e0/
c sq2pil = alog(sqrt(2.*pi)),  sqpi2l = alog (sqrt(pi/2.))
      data sqpi2l / 0.2257913526 4472743e0/
      data pi     / 3.1415926535 8979324e0/
c
      data xmax, dxrel / 0., 0. /
c
      if (xmax.ne.0.) go to 10
      xmax = r1mach(2)/alog(r1mach(2))
      dxrel = sqrt (r1mach(4))
c
 10   y = abs(x)
      if (y.gt.10.0) go to 20
c
c alog (abs (gamma(x))) for  abs(x) .le. 10.0
c
      alngam = alog (abs (gamma(x)))
      return
c
c alog (abs (gamma(x))) for abs(x) .gt. 10.0
c
 20   continue
c      if (y.gt.xmax) call seteru (
c     1  38halngam  abs(x) so big alngam overflows, 38, 2, 2)
c
      if (x.gt.0.) alngam = sq2pil + (x-0.5)*alog(x) - x + r9lgmc(y)
      if (x.gt.0.) return
c
      sinpiy = abs (sin(pi*y))
c      if (sinpiy.eq.0.) call seteru (31halngam  x is a negative integer,
c     1  31, 3, 2)
c
      alngam = sqpi2l + (x-0.5)*alog(y) - x - alog(sinpiy) - r9lgmc(y)
c
c      if (abs((x-aint(x-0.5))*alngam/x).lt.dxrel) call seteru (
c     1'alngam  answer lt half precision because x too near negative ',
c     2   68, 1, 1)
      return
c
      end
      function alnrel (x)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
      dimension alnrcs(23)
      external alog, csevl, inits, r1mach
c
c series for alnr       on the interval -3.75000d-01 to  3.75000d-01
c                                        with weighted error   1.93e-17
c                                         log weighted error  16.72
c                               significant figures required  16.44
c                                    decimal places required  17.40
c
      data alnrcs( 1) /   1.0378693562 743770e0 /
      data alnrcs( 2) /   -.1336430150 4908918e0 /
      data alnrcs( 3) /    .0194082491 35520563e0 /
      data alnrcs( 4) /   -.0030107551 12753577e0 /
      data alnrcs( 5) /    .0004869461 47971548e0 /
      data alnrcs( 6) /   -.0000810548 81893175e0 /
      data alnrcs( 7) /    .0000137788 47799559e0 /
      data alnrcs( 8) /   -.0000023802 21089435e0 /
      data alnrcs( 9) /    .0000004164 04162138e0 /
      data alnrcs(10) /   -.0000000735 95828378e0 /
      data alnrcs(11) /    .0000000131 17611876e0 /
      data alnrcs(12) /   -.0000000023 54670931e0 /
      data alnrcs(13) /    .0000000004 25227732e0 /
      data alnrcs(14) /   -.0000000000 77190894e0 /
      data alnrcs(15) /    .0000000000 14075746e0 /
      data alnrcs(16) /   -.0000000000 02576907e0 /
      data alnrcs(17) /    .0000000000 00473424e0 /
      data alnrcs(18) /   -.0000000000 00087249e0 /
      data alnrcs(19) /    .0000000000 00016124e0 /
      data alnrcs(20) /   -.0000000000 00002987e0 /
      data alnrcs(21) /    .0000000000 00000554e0 /
      data alnrcs(22) /   -.0000000000 00000103e0 /
      data alnrcs(23) /    .0000000000 00000019e0 /
c
      data nlnrel, xmin /0, 0./
c
      if (nlnrel.ne.0) go to 10
      nlnrel = inits (alnrcs, 23, 0.1*r1mach(3))
      xmin = -1.0 + sqrt(r1mach(4))
c
 10   continue
c      if (x.le.(-1.0)) call seteru (
c     1  18halnrel  x is le -1, 18, 2, 2)
c      if (x.lt.xmin) call seteru (
c     1  54halnrel  answer lt half precision because x too near -1, 54,
c     2  1, 1)
c
      if (abs(x).le.0.375) alnrel = x*(1. -
     1  x*csevl (x/.375, alnrcs, nlnrel))
      if (abs(x).gt.0.375) alnrel = alog (1.0+x)
c
      return
      end
      function gamma (x)
c jan 1984 edition.   w. fullerton, c3, los alamos scientific lab.
      dimension gcs(23)
      external alog, csevl, inits, r1mach, r9lgmc
c
      data gcs   ( 1) / .0085711955 90989331e0/
      data gcs   ( 2) / .0044153813 24841007e0/
      data gcs   ( 3) / .0568504368 1599363e0/
      data gcs   ( 4) /-.0042198353 96418561e0/
      data gcs   ( 5) / .0013268081 81212460e0/
      data gcs   ( 6) /-.0001893024 529798880e0/
      data gcs   ( 7) / .0000360692 532744124e0/
      data gcs   ( 8) /-.0000060567 619044608e0/
      data gcs   ( 9) / .0000010558 295463022e0/
      data gcs   (10) /-.0000001811 967365542e0/
      data gcs   (11) / .0000000311 772496471e0/
      data gcs   (12) /-.0000000053 542196390e0/
      data gcs   (13) / .0000000009 193275519e0/
      data gcs   (14) /-.0000000001 577941280e0/
      data gcs   (15) / .0000000000 270798062e0/
      data gcs   (16) /-.0000000000 046468186e0/
      data gcs   (17) / .0000000000 007973350e0/
      data gcs   (18) /-.0000000000 001368078e0/
      data gcs   (19) / .0000000000 000234731e0/
      data gcs   (20) /-.0000000000 000040274e0/
      data gcs   (21) / .0000000000 000006910e0/
      data gcs   (22) /-.0000000000 000001185e0/
      data gcs   (23) / .0000000000 000000203e0/
c
      data pi /3.14159 26535 89793 24e0/
c sq2pil is alog (sqrt (2.*pi) )
      data sq2pil /0.91893 85332 04672 74e0/
      data ngcs, xmin, xmax, xsml, dxrel /0, 4*0.0 /
c
      if (ngcs.ne.0) go to 10
c
c ---------------------------------------------------------------------
c initialize.  find legal bounds for x, and determine the number of
c terms in the series required to attain an accuracy ten times better
c than machine precision.
c
      ngcs = inits (gcs, 23, 0.1*r1mach(3))
c
      call r9gaml (xmin, xmax)
      xsml = exp (amax1 (alog (r1mach(1)), -alog(r1mach(2))) + 0.01)
      dxrel = sqrt (r1mach(4))
c
c ---------------------------------------------------------------------
c finish initialization.  start evaluating gamma(x).
c
 10   y = abs(x)
      if (y.gt.10.0) go to 50
c
c compute gamma(x) for abs(x) .le. 10.0.  reduce interval and
c find gamma(1+y) for 0. .le. y .lt. 1. first of all.
c
      n = x
      if (x.lt.0.) n = n - 1
      y = x - float(n)
      n = n - 1
      gamma = 0.9375 + csevl(2.*y-1., gcs, ngcs)
      if (n.eq.0) return
c
      if (n.gt.0) go to 30
c
c compute gamma(x) for x .lt. 1.
c
      n = -n
c      if (x.eq.0.) call seteru (14hgamma   x is 0, 14, 4, 2)
c      if (x.lt.0. .and. x+float(n-2).eq.0.) call seteru (
c     1  31hgamma   x is a negative integer, 31, 4, 2)
c      if (x.lt.(-0.5) .and. abs((x-aint(x-0.5))/x).lt.dxrel) call
c     1  seteru (68hgamma   answer lt half precision because x too near n
c     2egative integer, 68, 1, 1)
c      if (y.lt.xsml) call seteru (
c     1  54hgamma   x is so close to 0.0 that the result overflows,
c     2  54, 5, 2)
c
      do 20 i=1,n
        gamma = gamma / (x+float(i-1))
 20   continue
      return
c
c gamma(x) for x .ge. 2.
c
 30   do 40 i=1,n
        gamma = (y+float(i))*gamma
 40   continue
      return
c
c compute gamma(x) for abs(x) .gt. 10.0.  recall y = abs(x).
c
 50   continue
c      if (x.gt.xmax) call seteru (32hgamma   x so big gamma overflows,
c     1  32, 3, 2)
c
      gamma = 0.
c      if (x.lt.xmin) call seteru (35hgamma   x so small gamma underflows
c     1  , 35, 2, 0)
      if (x.lt.xmin) return
c
      gamma = exp((y-0.5)*alog(y) - y + sq2pil + r9lgmc(y) )
      if (x.gt.0.) return
c
c      if (abs((x-aint(x-0.5))/x).lt.dxrel) call seteru (
c     1  61hgamma   answer lt half precision, x too near negative integer
c     2  , 61, 1, 1)
c
      sinpiy = sin (pi*y)
c      if (sinpiy.eq.0.) call seteru (
c     1  31hgamma   x is a negative integer, 31, 4, 2)
c
      gamma = -pi / (y*sinpiy*gamma)
c
      return
      end
      function r9lgmc (x)
c august 1977 edition.  w. fullerton, c3, los alamos scientific lab.
c
c compute the log gamma correction factor for x .ge. 10.0 so that
c  alog (gamma(x)) = alog(sqrt(2*pi)) + (x-.5)*alog(x) - x + r9lgmc(x)
c
      dimension algmcs(6)
      external alog, csevl, inits, r1mach
c
c series for algm       on the interval  0.          to  1.00000d-02
c                                        with weighted error   3.40e-16
c                                         log weighted error  15.47
c                               significant figures required  14.39
c                                    decimal places required  15.86
c
      data algmcs( 1) /    .1666389480 45186e0 /
      data algmcs( 2) /   -.0000138494 817606e0 /
      data algmcs( 3) /    .0000000098 108256e0 /
      data algmcs( 4) /   -.0000000000 180912e0 /
      data algmcs( 5) /    .0000000000 000622e0 /
      data algmcs( 6) /   -.0000000000 000003e0 /
c
      data nalgm, xbig, xmax / 0, 2*0.0 /
c
      if (nalgm.ne.0) go to 10
      nalgm = inits (algmcs, 6, r1mach(3))
      xbig = 1.0/sqrt(r1mach(3))
      xmax = exp (amin1(alog(r1mach(2)/12.0), -alog(12.0*r1mach(1))) )
c
 10   continue
c      if (x.lt.10.0) call seteru (23hr9lgmc  x must be ge 10, 23, 1, 2)
      if (x.ge.xmax) go to 20
c
      r9lgmc = 1.0/(12.0*x)
      if (x.lt.xbig) r9lgmc = csevl (2.0*(10./x)**2-1., algmcs, nalgm)/x
      return
c
 20   r9lgmc = 0.0
c      call seteru (34hr9lgmc  x so big r9lgmc underflows, 34, 2, 0)
      return
c
      end
      function csevl (x, cs, n)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c evaluate the n-term chebyshev series cs at x.  adapted from
c r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973).  also see fox
c and parker, chebyshev polys in numerical analysis, oxford press, p.56.
c
c             input arguments --
c x      value at which the series is to be evaluated.
c cs     array of n terms of a chebyshev series.  in eval-
c        uating cs, only half the first coef is summed.
c n      number of terms in array cs.
c
      dimension cs(n)
c
c      if (n.lt.1) call seteru (28hcsevl   number of terms le 0, 28, 2,2)
c      if (n.gt.1000) call seteru (31hcsevl   number of terms gt 1000,
c     1  31, 3, 2)
c      if (x.lt.(-1.1) .or. x.gt.1.1) call seteru (
c     1  25hcsevl   x outside (-1,+1), 25, 1, 1)
c
      b1 = 0.
      b0 = 0.
      twox = 2.*x
      do 10 i=1,n
        b2 = b1
        b1 = b0
        ni = n + 1 - i
        b0 = twox*b1 - b2 + cs(ni)
 10   continue
c
      csevl = 0.5 * (b0-b2)
c
      return
      end
      function inits (os, nos, eta)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c initialize the orthogonal series so that inits is the number of terms
c needed to insure the error is no larger than eta.  ordinarily, eta
c will be chosen to be one-tenth machine precision.
c
c             input arguments --
c os     array of nos coefficients in an orthogonal series.
c nos    number of coefficients in os.
c eta    requested accuracy of series.
c
      dimension os(nos)
c
c      if (nos.lt.1) call seteru (
c     1  35hinits   number of coefficients lt 1, 35, 2, 2)
c
      err = 0.
      do 10 ii=1,nos
        i = nos + 1 - ii
        err = err + abs(os(i))
        if (err.gt.eta) go to 20
 10   continue
c
 20   continue
c      if (i.eq.nos) call seteru (28hinits   eta may be too small, 28,
c     1  1, 2)
      inits = i
c
      return
      end
      subroutine r9gaml (xmin, xmax)
c april 1977 version.  w. fullerton, c3, los alamos scientific lab.
c
c calculate the minimum and maximum legal bounds for x in gamma(x).
c xmin and xmax are not the only bounds, but they are the only non-
c trivial ones to calculate.
c
c             output arguments --
c xmin   minimum legal value of x in gamma(x).  any smaller value of
c        x might result in underflow.
c xmax   maximum legal value of x in gamma(x).  any larger value will
c        cause overflow.
c
      external alog, r1mach
c
      alnsml = alog(r1mach(1))
      xmin = -alnsml
      do 10 i=1,10
        xold = xmin
        xln = alog(xmin)
        xmin = xmin - xmin*((xmin+0.5)*xln - xmin - 0.2258 + alnsml)
     1    / (xmin*xln + 0.5)
        if (abs(xmin-xold).lt.0.005) go to 20
 10   continue
c      call seteru (27hr9gaml  unable to find xmin, 27, 1, 2)
c
 20   xmin = -xmin + 0.01
c
      alnbig = alog(r1mach(2))
      xmax = alnbig
      do 30 i=1,10
        xold = xmax
        xln = alog(xmax)
        xmax = xmax - xmax*((xmax-0.5)*xln - xmax + 0.9189 - alnbig)
     1    / (xmax*xln - 0.5)
        if (abs(xmax-xold).lt.0.005) go to 40
 30   continue
c      call seteru (27hr9gaml  unable to find xmax, 27, 2, 2)
c
 40   xmax = xmax - 0.01
      xmin = amax1 (xmin, -xmax+1.)
c
      return
      end
      subroutine r9upak (x, y, n)
c august 1980 portable edition.  w. fullerton, los alamos scientific lab
c
c unpack floating point number x so that x = y * 2.0**n, where
c 0.5 .le. abs(y) .lt. 1.0 .
c
      absx = abs(x)
      n = 0
      y = 0.0
      if (x.eq.0.0) return
c
 10   if (absx.ge.0.5) go to 20
      n = n - 1
      absx = absx*2.0
      go to 10
c
 20   if (absx.lt.1.0) go to 30
      n = n + 1
      absx = absx*0.5
      go to 20
c
 30   y = sign (absx, x)
      return
c
      end
      subroutine tautsp(tau,gtau,ntau,gamma,s,break,coef,l,k,iflag)
c  constructs cubic spline interpolant to given data
c       tau(i),gtau(i),i=1,...ntau
c  if gamma.gt.0 additional knots are introduced where needed to
c  make the interpolant more flexible locally. this avoids extraneous
c  inflection points typical of cubic spline interpolation at knots to
c  rapidly changing data.
c
c  taken from c.de boor, a practical guide to splines, pp 310-314
c  modified to remove error print statement, a.chave igpp mar 1982
c
c  parameters
c
c      input
c  tau     sequence of data points, monotonically increasing
c  gtau    corresponding sequence of function values
c  ntau    number of data points, must be .ge.4
c  gamma   indicates whether additional flexibility is required
c          =0., no additional knots
c          =(0.,3.), under certain conditions on the given data at
c                  points i-1,...i+2 a knot is added in the i-th
c                  interval, i=2,...ntau-2. see description below.
c                  the interpolant gets rounded with increasing gamma.
c                  a value of 2.5 is typical.
c          =(3.,6.), same, except that knots might also be added in
c                  intervals in which an inflection point would be
c                  permitted. a value of 5.5 is typical.
c      output
c  break,coef,l,k give the pp-representation of the interpolant
c      for break(i).le.x.le.break(i+1) the interpolant has the form
c      f(x)=coef(1,i)+dx(coef(2,i)+(dx/2)(coef(3,i)+(dx/3)coef(4,i)))
c      with dx=x-break(i) and i=1,...l
c  iflag=0, ok
c       =1, input was incorrect
c      workspace
c  s      is of size (ntau,6). the individual columns of this array
c         contain
c         s(.,1)=dtau=tau(.+1)-tau(.)
c         s(.,2)=diag=diagonal in linear system
c         s(.,3)=u=upper diagonal in linear system
c         s(.,4)=r=right side for linear system(initially)
c         s(.,5)=z=indicator of additional knots
c         s(.,6)=1/hsecnd(1,x) with x=z or 1-z, see below
c
c  *****method*****
c  on the i-th interval, (tau(i),tau(i+1)), the interpolant is of the
c  form
c    f(u(x))=a+b*u+c*h(u,z)+d*h(1-u,1-z)
c  with u(x)=(x-tau(i))/dtau(i), where
c    z(i)=addg(i+1)/(addg(i)+addg(i+1))
c    (=.5 in the case where the denominator vanishes)
c    addg(j)=abs(ddg(j)),  ddg(j)=dg(j+1)-dg(j),
c    dg(j)=divdif(j)=(gtau(j+1)-gtau(j))/dtau(j)
c  and
c    h(u,z)=alpha*u**3+(1-alpha)*(max(((u-zeta)/(1-zeta)),0)**3
c  with
c    alpha(z)=(1-gamma/3)/zeta
c    zeta(z)=1-gamma*min((1-z),1/3)
c  for 1/3.le.z.le.2/3 f is a cubic polynomial on the interval i.
c  otherwise it has one additional knot at
c    tau(i)+zeta*dtau(i)
c  as z approaches 1, h(.,z) has an increasingly sharp bend near 1,
c  allowing f to turn rapidly near the additional knot.
c  in terms of f(j)=gtau(j) and fsecnd(j)=2 derivative of f at tau(j),
c  the coefficients for f(u(x)) are given as
c    a=f(i)-d
c    b=(f(i+1)-f(i))-(c-d)
c    c=fsecnd(i+1)*dtau(i)**2/hsecnd(1,z)
c    d=fsecnd(i)*dtau(i)**2/hsecnd(1,1-z)
c  hence can be computed once fsecnd(i),i=1,...ntau is fixed.
c  f is automatically continuous and has a continuous second
c  derivative, except when z=0 or 1 for some i. we determine
c  fsecnd(.) from the requirement that the first derivative of
c  f be continuous. in addition we require that the third derivative
c  be continuous across tau(2) and across tau(ntau-1). this leads to
c  a strictly diagonally dominant tridiagonal linear system for the
c  fsecnd(i) which we solve by gauss elimination without pivoting.
c
      real break(*),coef(4,*),gamma,gtau(*),s(ntau,6),tau(*),
     $       alpha,c,d,del,denom,divdif,entry,entry3,factor,factr2,gam,
     $       onemzt,ratio,sixth,temp,x,z,zeta,zt2
      integer n
      alph(x)=min(1.,onemg3/x)
c  there must be at least 4 interpolation points
      if(ntau.lt.4)then
        iflag=1
        return
      endif
      ntaum1=ntau-1
      do 6 i=1,ntaum1
        s(i,1)=tau(i+1)-tau(i)
        if(s(i,1).le.0.)then
c  data out of order-error return
          iflag=1
          return
        endif
   6    s(i+1,4)=(gtau(i+1)-gtau(i))/s(i,1)
      do 7 i=2,ntaum1
   7    s(i,4)=s(i+1,4)-s(i,4)
c
c  construct system of equations for second derivatives at tau. at each
c  interior data point there is one continuity equation, at the first
c  and the last interior data point there is an additional one for a
c  total of ntau equations in ntau unknowns.
c
      i=2
      s(2,2)=s(1,1)/3.
      sixth=1./6.
      method=2
      gam=gamma
      if(gam.le.0.)method=1
      if(gam.le.3.)goto 9
      method=3
      gam=gam-3.
    9 onemg3=1.-gam/3.
c  loop over i
   10 continue
c  construct z(i) and zeta(i)
      z=.5
      goto (19,11,12)method
   11 if(s(i,4)*s(i+1,4).lt.0.)goto 19
   12 temp=abs(s(i+1,4))
      denom=abs(s(i,4))+temp
      if(denom.eq.0.)goto 19
      z=temp/denom
      if(abs(z-.5).le.sixth)z=.5
   19 s(i,5)=z
c  set up the part of the i-th equation which depends on the
c  i-th interval
      if(z-.5)21,22,23
   21 zeta=gam*z
      onemzt=1.-zeta
      zt2=zeta*zeta
      alpha=alph(onemzt)
      factor=zeta/(alpha*(zt2-1.)+1.)
      s(i,6)=zeta*factor/6.
      s(i,2)=s(i,2)+s(i,1)*((1.-alpha*onemzt)*factor/2.-
     $       s(i,6))
c  if z=0 and the previous z=1, then d(i)=0, since then
c  u(i-1)=l(i+1)=0 and its value doesn't matter. reset d(i)=1
c  to insure nonzero pivot in elimination.
      if(s(i,2).le.0.)s(i,2)=1.
      s(i,3)=s(i,1)/6.
      goto 25
   22 s(i,2)=s(i,2)+s(i,1)/3.
      s(i,3)=s(i,1)/6.
      goto 25
   23 onemzt=gam*(1.-z)
      zeta=1.-onemzt
      alpha=alph(zeta)
      factor=onemzt/(1.-alpha*zeta*(1.+onemzt))
      s(i,6)=onemzt*factor/6.
      s(i,2)=s(i,2)+s(i,1)/3.
      s(i,3)=s(i,6)*s(i,1)
   25 if(i.gt.2)goto 30
      s(1,5)=.5
c  the first two equations enforce continuity of the first and third
c  derivatives across tau(2)
      s(1,2)=s(1,1)/6.
      s(1,3)=s(2,2)
      entry3=s(2,3)
      if(z-.5)26,27,28
   26 factr2=zeta*(alpha*(zt2-1.)+1.)/(alpha*(zeta*zt2-1.)+1.)
      ratio=factr2*s(2,1)/s(1,2)
      s(2,2)=factr2*s(2,1)+s(1,1)
      s(2,3)=-factr2*s(1,1)
      goto 29
   27 ratio=s(2,1)/s(1,2)
      s(2,2)=s(2,1)+s(1,1)
      s(2,3)=-s(1,1)
      goto 29
   28 ratio=s(2,1)/s(1,2)
      s(2,2)=s(2,1)+s(1,1)
      s(2,3)=-s(1,1)*6.*alpha*s(2,6)
c  at this point the first two equations read
c    diag(1)*x1+u(1)*x2+entry3*x3=r(2)
c    -ratio*diag(1)*x1+diag(2)*x2+u(2)*x3=0
c  eliminate first unknown from second equation
   29 s(2,2)=ratio*s(1,3)+s(2,2)
      s(2,3)=ratio*entry3+s(2,3)
      s(1,4)=s(2,4)
      s(2,4)=ratio*s(1,4)
      goto 35
   30 continue
c  the i-th equation enforces continuity of the first derivative across
c  tau(i). it has been set up in statements 35 to 40 and 21 to 25 and
c  reads
c     -ratio*diag(i-1)*xi-1+diag(i)*xi+u(i)*xi+1=r(i)
c  eliminate (i-1)st unknown from this equation
      s(i,2)=ratio*s(i-1,3)+s(i,2)
      s(i,4)=ratio*s(i-1,4)+s(i,4)
c
c  set up the part of the next equation which depends on the i-th
c  interval
   35 if(z-.5)36,37,38
   36 ratio=-s(i,6)*s(i,1)/s(i,2)
      s(i+1,2)=s(i,1)/3.
      goto 40
   37 ratio=-(s(i,1)/6.)/s(i,2)
      s(i+1,2)=s(i,1)/3.
      goto 40
   38 ratio=-(s(i,1)/6.)/s(i,2)
      s(i+1,2)=s(i,1)*((1.-zeta*alpha)*factor/2.-s(i,6))
c end of i loop
   40 i=i+1
      if(i.lt.ntaum1)goto 10
      s(i,5)=.5
c
c  last two equations
c  the last two equations enforce continuity of the first and third
c  derivatives across tau(ntau-1)
      entry=ratio*s(i-1,3)+s(i,2)+s(i,1)/3.
      s(i+1,2)=s(i,1)/6.
      s(i+1,4)=ratio*s(i-1,4)+s(i,4)
      if(z-.5)41,42,43
   41 ratio=s(i,1)*6.*s(i-1,6)*alpha/s(i-1,2)
      s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
      s(i,3)=-s(i-1,1)
      goto 45
   42 ratio=s(i,1)/s(i-1,2)
      s(i,2)=ratio*s(i-1,3)+s(i,1)+s(i-1,1)
      s(i,3)=-s(i-1,1)
      goto 45
   43 factr2=onemzt*(alpha*(onemzt**2-1.)+1.)/
     $       (alpha*(onemzt**3-1.)+1.)
      ratio=factr2*s(i,1)/s(i-1,2)
      s(i,2)=ratio*s(i-1,3)+factr2*s(i-1,1)+s(i,1)
      s(i,3)=-factr2*s(i-1,1)
c  at this point the last two equations read
c    diag(i)*xi+u(i)*xi+1=r(i)
c    -ratio*diag(i)*xi+diag(i+1)*xi+1=r(i)
c  eliminate xi from the last equation
   45 s(i,4)=ratio*s(i-1,4)
      ratio=-entry/s(i,2)
      s(i+1,2)=ratio*s(i,3)+s(i+1,2)
      s(i+1,4)=ratio*s(i,4)+s(i+1,4)
c
c back substitution
c
      s(ntau,4)=s(ntau,4)/s(ntau,2)
   50 s(i,4)=(s(i,4)-s(i,3)*s(i+1,4))/s(i,2)
      i=i-1
      if(i.gt.1)goto 50
      s(1,4)=(s(1,4)-s(1,3)*s(2,4)-entry3*s(3,4))/s(1,2)
c
c  construct polynomial pieces
c
      break(1)=tau(1)
      l=1
      do 70 i=1,ntaum1
        coef(1,l)=gtau(i)
        coef(3,l)=s(i,4)
        divdif=(gtau(i+1)-gtau(i))/s(i,1)
        z=s(i,5)
        if(z-.5)61,62,63
   61   if(z.eq.0.)goto 65
        zeta=gam*z
        onemzt=1.-zeta
        c=s(i+1,4)/6.
        d=s(i,4)*s(i,6)
        l=l+1
        del=zeta*s(i,1)
        break(l)=tau(i)+del
        zt2=zeta*zeta
        alpha=alph(onemzt)
        factor=onemzt**2*alpha
        coef(1,l)=gtau(i)+divdif*del+s(i,1)**2*(d*onemzt*
     $            (factor-1.)+c*zeta*(zt2-1.))
        coef(2,l)=divdif+s(i,1)*(d*(1.-3.*factor)+c*(3.*
     $            zt2-1.))
        coef(3,l)=6.*(d*alpha*onemzt+c*zeta)
        coef(4,l)=6.*(c-d*alpha)/s(i,1)
        coef(4,l-1)=coef(4,l)-6.*d*(1.-alpha)/(del*zt2)
        coef(2,l-1)=coef(2,l)-del*(coef(3,l)-del/2.*
     $              coef(4,l-1))
        goto 68
   62   coef(2,l)=divdif-s(i,1)*(2.*s(i,4)+s(i+1,4))/6.
        coef(4,l)=(s(i+1,4)-s(i,4))/s(i,1)
        goto 68
   63   onemzt=gam*(1.-z)
        if(onemzt.eq.0.)goto 65
        zeta=1.-onemzt
        alpha=alph(zeta)
        c=s(i+1,4)*s(i,6)
        d=s(i,4)/6.
        del=zeta*s(i,1)
        break(l+1)=tau(i)+del
        coef(2,l)=divdif-s(i,1)*(2.*d+c)
        coef(4,l)=6.*(c*alpha-d)/s(i,1)
        l=l+1
        coef(4,l)=coef(4,l-1)+6.*(1.-alpha)*c/(s(i,1)*onemzt**3)
        coef(3,l)=coef(3,l-1)+del*coef(4,l-1)
        coef(2,l)=coef(2,l-1)+del*(coef(3,l-1)+del/2.*coef(4,l-1))
        coef(1,l)=coef(1,l-1)+del*(coef(2,l-1)+del/2.*(coef(3,l-1)
     $            +del/3.*coef(4,l-1)))
        goto 68
   65   coef(2,l)=divdif
        coef(3,l)=0.
        coef(4,l)=0.
   68   l=l+1
   70   break(l)=tau(i+1)
      l=l-1
      k=4
      iflag=0
      return
      end
      real function ppvalu(break,coef,l,k,x,jderiv)
c calculates value at x of the jderiv-th derivative of pp fct
c from pp-representation
c
c  break,coef,l,k are as output by tautsp
c  x is the value of x at which the function is to be evaluated
c  jderiv is an integer .ge.0 giving the derviative to be evaluated
      real break(l),coef(k,l),x,fmmjdr,h
      ppvalu=0.
      fmmjdr=k-jderiv
c  derivatives of order .ge. k are zero
      if(fmmjdr.le.0.)return
c  find index i of largest breakpoint to the left of x
      call interv(break,l,x,i,ndummy)
c  evaluate jderiv-th derivative of i-th polynomial piece at x
      h=x-break(i)
      do 10 m=k,jderiv+1,-1
        ppvalu=(ppvalu/fmmjdr)*h+coef(m,i)
   10   fmmjdr=fmmjdr-1.
      return
      end
      subroutine interv(xt,lxt,x,left,mflag)
c computes left=max(i,1.le.i.le..lxt .and. xt(i).le.x)
      real x,xt(*)
      data ilo/1/
      save ilo
      ihi=ilo+1
      if(ihi.lt.lxt)goto 20
      if(x.ge.xt(lxt))goto 110
      if(lxt.le.1)goto 90
      ilo=lxt-1
      ihi=lxt
   20 if(x.ge.xt(ihi))goto 40
      if(x.ge.xt(ilo))goto 100
      istep=1
   31 ihi=ilo
      ilo=ihi-istep
      if(ilo.le.1)goto 35
      if(x.ge.xt(ilo))goto 50
      istep=istep*2
      goto 31
   35 ilo=1
      if(x.lt.xt(1))goto 90
      goto 50
   40 istep=1
   41 ilo=ihi
      ihi=ilo+istep
      if(ihi.ge.lxt)goto 45
      if(x.lt.xt(ihi))goto 50
      istep=istep*2
      goto 41
   45 if(x.ge.xt(lxt))goto 110
      ihi=lxt
   50 middle=(ilo+ihi)/2
      if(middle.eq.ilo)goto 100
      if(x.lt.xt(middle))goto 53
      ilo=middle
      goto 50
   53 ihi=middle
      goto 50
   90 mflag=-1
      left=1
      return
  100 mflag=0
      left=ilo
      return
  110 mflag=1
      left=lxt
      return
      end

      real function xmedian(x,n)
      double precision sum
      dimension x(n)
      external fmedian
      sum=0.d0
      do 10 i=1,n
        sum=sum+x(i)
        b=min(b,x(i))
        c=max(c,x(i))
  10    continue
      r=sum/n
      call fzero(fmedian,b,c,r,1.e-6,1.e-32,iflag,x,[n])
      xmedian=b
      return
      end

      real function fmedian(xmed,x,n)
      double precision sum
      dimension x(n)
      sum=0.d0
      do 10 i=1,n
        sum=sum+(x(i)-xmed)/abs(x(i)-xmed)
  10    continue
      fmedian=sum
      return
      end
      
      
