Solved: Certainly not - I was just

milenko1976 · ‎09-14-2016

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?

mecej4 · ‎09-15-2016

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

mecej4 · ‎09-14-2016

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?

TimP · ‎09-14-2016

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.

Steven_L_Intel1 · ‎09-14-2016

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

TimP · ‎09-14-2016

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).

Steven_L_Intel1 · ‎09-14-2016

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

milenko1976 · ‎09-15-2016

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.

mecej4 · ‎09-15-2016

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?

milenko1976 · ‎09-15-2016

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

milenko1976 · ‎09-15-2016

I am attaching the whole folder

524077

mecej4 · ‎09-15-2016

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.

mecej4 · ‎09-15-2016

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.

milenko1976 · ‎09-16-2016

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.