Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Highlighted
Beginner
24 Views

If the actual argument is scalar, the dummy argument shall be scalar

Jump to solution

Following Steve Lionel'sadvicee with warn interface I got this warning for the code I am compiling

math.f(2079): error #8284: If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic.   [IPAR]

 

In my subroutine fzero I have vector if integer parameters ipar.The guy who wrote the code declared

 dimension rpar(*),ipar(*)

 

How should I change this to function properly,not just to silence the compiler?

 

0 Kudos

Accepted Solutions
Highlighted
Black Belt
24 Views

Here is a simple fix. In the

Jump to solution

Here is a simple fix. In the caller of FZERO, declare an integer array nn(1), set nn(1) = n and pass nn as the last argument to FZERO. In this way, you will satisfy the formal requirement that no scalar can be passed as the actual argument when the corresponding formal argument is an array.

 

View solution in original post

0 Kudos
12 Replies
Highlighted
Black Belt
24 Views

We need to be shown the call

Jump to solution

We need to be shown the call statement as well as the subroutine heading, along with the declarations of all the variables used in the call.

Is the argument IPAR actually passed and used, or is this a place-holder argument?

0 Kudos
Highlighted
Black Belt
24 Views

If those arguments are

Jump to solution

If those arguments are scalars in the caller (in every use), you should remove the (*) in the callee.  If both scalars and arrays are used in the callers, the scalar can be made into an array (length 1) by placing [] around it in the CALL.  The difference is flagged by ifort only as a standards compliance helper; this violation of standard "works" in most compilers for IA targets (if there is no difference in data types).

Fortran standard doesn't allow assumed size dummy arrays to be compatible with scalars, except possibly as an extension, which the compiler must be capable of flagging.

I'd have some concerns about allowing the compiler to implement inter-procedural optimization with some of the wilder cases you have put forward.  You might try to make it work with -O0 and with -fno-inline-functions before fixing those problems.

0 Kudos
Highlighted
24 Views

See the section on Sequence

Jump to solution

See the section on Sequence Association in https://software.intel.com/en-us/blogs/2009/03/31/doctor-fortran-in-ive-come-here-for-an-argument

Retired 12/31/2016
0 Kudos
Highlighted
Black Belt
24 Views

I suppose sequence

Jump to solution

I suppose sequence association must have been valid for the application to work at all in the past, but I don't think it can be trusted for cases such as mixing real((8) and integer(4).

0 Kudos
Highlighted
24 Views

Certainly not - I was just

Jump to solution

Certainly not - I was just addressing the specific issue in this thread.

Retired 12/31/2016
0 Kudos
Highlighted
Beginner
24 Views

My program is here

Jump to solution

My program is here

https://app.box.com/files/0/f/0/1/f_94899963580

 

I have tried to compile like this

FCFLAGS= -O0 -fno-inline-functions -warn interface

but everything is the same.Tim p. has mentioned that I should remove(*).I have not done that because the guy who wrote the code was having in mind assumed value array at least,as I understand.

0 Kudos
Highlighted
Black Belt
24 Views

The Box link appears to be

Jump to solution

The Box link appears to be for your private files and the target does not have public access. Why not simply attach the source file here?

0 Kudos
Highlighted
Beginner
24 Views

math.f

      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(*),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(*)
      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
      
      

0 Kudos
Highlighted
Beginner
24 Views

I am attaching the whole

Jump to solution

I am attaching the whole folder

524077

0 Kudos
Highlighted
Black Belt
25 Views

Here is a simple fix. In the

Jump to solution

Here is a simple fix. In the caller of FZERO, declare an integer array nn(1), set nn(1) = n and pass nn as the last argument to FZERO. In this way, you will satisfy the formal requirement that no scalar can be passed as the actual argument when the corresponding formal argument is an array.

 

View solution in original post

0 Kudos
Highlighted
Black Belt
24 Views

When I looked at your code in

Jump to solution

When I looked at your code in order to answer your Fortran-related question, I noticed some bugs that may be worthwhile to investigate.

  • In function XMEDIAN, the variables b and c are undefined when first used in the DO loop.
  • In function FMEDIAN, you compute (x(i)-xmed)/abs(x(i)-xmed). If the array x contains an odd number of elements, the median will be, by definition, one of the elements of the array, and the expression will be evaluated as 0/|0|. Apart from this problem, note that the expression will have only the values +1 or -1 if the candidate xmed is not a member of the set x(1..n). Thus, you could compute fmedian as count(x > xmed) - count(x < xmed).
  • The code finds the median of an array x(1..n) by calling the nonlinear zero finding routine FZERO from the SLATEC library. This is probably a poor and incorrect method for several reasons: (a) FZERO is for finding the zero of a continuous function of a real variable. Your function is a discrete function, and (b) the solution returned by using FZERO does not satisfy the definition of the median of an array of values.
  • Look for the Blum-Floyd-Rivest-Tarjan algorithm or the slightly simpler QuickSelect algorithm if you want to find the true median.
0 Kudos
Highlighted
Beginner
24 Views

Thanks mecej4,The code is not

Jump to solution

Thanks mecej4,The code is not mine and I am trying just to compile it to fit for the purpose.Yes there are some bugs,I agree.

0 Kudos