INDEX of routines in the accumulating mathematical library (January 2006)

BS 	completes the L2 soln. of overdetermined linear equations after HH
CFGEN	generates the 3-term recurrence cfts. for a summation inner producr
CFDIRZ	direct evaluation of a continued fraction with complex argument
DENIN	evaluates the integrated density of states from the 3-term rec.reln. 
DENQD	evaluates the density of states by the quadrature method
DOSCAL	computes the density of states from the 3-term rec.reln. 2005 method
DYCFTL	estimates decay rate of  Christoffel function after DOSCAL,GRCALZ
GRASSZ	estimates the Greenian of a J-fraction for complex arg asymptotic limit
GRCALZ	estimates the Greenian of a J-fraction for complex arg max current B
GRGENZ	estimates the Greenian of a J-fraction for complex arg max current A
HH	factorises a rectangular matrix prior to BS completing the L2 solution
INFGEN	generates a quadrature for integrals over the doubly infinite line
INITAL	tabulates quadrature points prior to use of the routine QFINT family
NUMC	counts no. of eigenvalues of a symm.3-diagonal matrix greater than x
NUMD	evaluates the inverse log derivative a the det of a tridiagonal matrix
PREPAS	prepares parameters prior to use of DYCFTL
QFINT	evaluates the definite integral of a regular function in a finite range
QVINT	evaluates a line integral of a regular function
QVINT1	evaluates the integral over a parallelogram of a regular function
QVINT2	evaluates the integral in a parallelepiped of a regular function
RECQD	calculates the cfts of the quadrature corresponding to a given J-matrix
RECRTS	computes some or all of the eigenvalues of a J-matrix
RECSTD	generates the inner product for some standard weight functions
RECWT	computes the weight at a node in a quadrature: Christoffel function
SIMBIS	finds a zero of a function by simple bisection
VECCF	evaluates a vector continued fraction
VECPOL	evaluates the vector polyns. And Christoffel fn. For a vector c.f



























      SUBROUTINE BS(A,P,MM,M,N,PT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MM,N),P(M),PT(N)
 
Finds the solution vector of the overdetermined system of linear
equations  Ax = p  by the method of least squares, after A has
been factorised by HH.
 
ARGUMENTS : (* denotes an overwritten argument)
 
     A	unaltered output from HH
     P*	on input, the right-hand vector of the equation
	on output, the solution vector is held in the first N places :
	the sum of squares of the last M-N components gives the sum of
	squares of deviations in the least squares problem
     MM	first DIMENSION of array A in calling program
     M	number of rows in matrix A
     N 	number of columns in matrix A
     PT	unaltered output from HH












      SUBROUTINE CFDIRZ(Z,A,B2,LEN,TERM,VALZ)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LEN),B2(LEN)
      COMPLEX*16 Z,TERM,VALZ

Evaluates a (J-) continued fraction with real coefficients and complex argument by direct use of the recurrence.

ARGUMENTS: (* denotes an overwritten argument)

   Z	complex argument of the J-fraction
   A	diagonal coefficients of the J-fraction
   B2	off-diagonal coefficients of the J-fraction
   LEN	number of levels in the J-fraction
   TERM	terminator value for the bottom of the J-fraction
   VALZ*	complex value of the J-fraction









      SUBROUTINE CFGEN(X,W,N,EPS,A,B2,LM1,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION X(N),W(N),A(LM1),B2(LM1),WK(N,2)
 
     Calculates the three-term recurrence relation corresponding to a weight function defined in terms of a summation inner product:

      <f,g> = sum i=1 to N { f(X(i)) * g( X(i)) * W(i) } .
 
Then, with p(x,0) = 1 and p(x,-1) = 0 , for n = 1,2...,LM1
 
      A(n) = <x * p(x,n-1) , p(x,n-1)> / < p(x,n-1) , p(x,n-1) >
 
     B2(n) = < p(x,n-1) , p(x,n-1) > / < p(x,n-2) , p(x,n-2) >
 
    p(x,n) = ( x-A(n) ) * p(x,n-1) - B2(n) * p(x,n-2)
 
The polynomials p(x,n) are represented by their values at the points X(i) , with an arbitrary normalisation to preserve numerical stability. [Haydock & Nex]
 
ARGUMENTS : (* indicates an overwritten argument)
 
     X	nodes of the summation inner product
     W	weights of the summation inner product
     N	number of terms in the summation inner product
     EPS	a value of B2(i) less than EPS terminates routine
     A*	diagonal coefficients in the recurrence A(i),i=1,LM1
     B2*	off-diagonal coefficients in the recurrence B2(i),i=1,LM1
     LM1*	on input  required length of recurrence 
	on output length of recurrence actually 
	generated
     WK	work array of length at least 2*N
























      FUNCTION DENIN(E,A,B2,NA,NP,LL,ALP,EPS,WK,IWK,ICODE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(NA,NP),B2(NA,NP),WK(NA,2),IWK(NA,2)
 
     Evaluates the integrated density of states, N(E), corresponding to a given sum of continued fractions (J-matrices)  at a given point E and returns that value, using the 'quadrature' approach.
This routine uses RECWT,RECRTS,NUMC,NUMD
 
Arguments : (* indicates an overwritten argument)
 
DENIN	takes the computed value of the integrated density of states at E
E	value at which integrated density of states required
A*	diagonal J-matrix elements (A(LL,k) overwritten) i=1,LL-1
B2	squares of sub-diagonal J-matrix elements i=2,LL
	B2(1,k) is the weight in the k th band
NA	first DIMENSION of arrays A , B2  , WK and IW >= LL
NP	number of density of states to be summed
LL	length of tridiagonalisations
ALP	proportion of weight at last node, 0<ALP<1 ,usually =0.5
EPS	accuracy required in root-finding
WK*	work array of length at least 2*NA
	first column contains nodes and second weights for the last
	 set of continued fraction coefficients
IWK*	integer array of length at least 2*NA (o/p from last RECRTS call)
ICODE*	on a successful output contains no.of nodes in last quadrature
	negative  if a multiple root in RECRTS (absolute value as above)





























      FUNCTION DENQD(E,EMX,A,B2,LL,ALP,EPS,TB,NT,NQ,NE,IWK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LL),B2(LL),TB(NT,4),IWK(LL)
 
     Evaluates the density of states, n(E), corresponding to a given continued fraction (J-matrix)  at a given point E and returns that value and also quadrature nodes and weights at a set of points bounded above by EMX. The table of values TB is output so that the integrated density of states, density of states , and similar functions may evaluated at each e(i) not greater than EMX. e.g. The integral to TB(i,1) of f(x).n(x) dx is approximated by
the sum j=1,i (last term times alpha) f(TB(j,1)) * TB(j,2).
The expressions defining the approximation are as follows (with n=LL):
 
A(n)= E - B2(n) * p(E,n-1)/p(E,n)
 
w(i) = q(e(i),n) / p'(e(i),n+1)
 
dw(i)/dA(LL) = the following expression evaluated at e(i)  :
 
{[q'(n)*p(n)-p'(n+1)*q(n-1)]+w(i)*[p'(n)*p'(n+1)-p"(n+1)*p(n)]} /
                                                          p'(n+1)**2
 
n(E) = p(E,n+1)/p(E,n)*{ [sum i<NE dw(i)/dA(LL)] + ALP* dw(NE)/dA(LL) }
 
p and q are the monic orthogonal polynomials of the first and
second kind associated with the weight function n(e) (see PLYVAL
for explicit definition of p ), and the e(i) are the eigenvalues
of the given Jacobi matrix with A(LL) appended so that e(NE)=E.
     In the actual calculation the values of the polynomials are
renormalsied to maintain numerical stability (only ratios of
polynomials appear in the above expressions).[Nex]
This routine uses RECWT,RECRTS,NUMC,NUMD
 
Arguments : (* indicates an overwritten argument)
 
DENQD	takes the computed value of the density of states at E
E	alue at which density of states required
EMX	upper limit of range of quadrature nodal values required > E
A*	diagonal J-matrix elements (A(LL) overwritten) i=1,LL-1
B2	squares of sub-diagonal J-matrix elements i=2,LL
B2(1)	is the total weight of the density of states
LL	length of tridiagonalisation
ALP	proportion of weight at last node, 0<ALP<1 ,usually =0.5
EPS	accuracy required in root-finding
TB*	table of quadrature nodes and differentials :
	TB(i,1)   nodal points : e(i)
	TB(i,2)   nodal weights : w(i)
	TB(i,3)   dw(i)/dA(LL)
	TB(i,4)   p'(e(i),LL+1) / p(e(i),LL)
NT	first DIMENSION of array TB in calling segment
NQ*	number of nodal values calculated
NE*	IABS(NE) gives index of node corresponding to E :TB(NE,1)=E
	if negative the accuracy is inadequate
	if =0 a multiple root was identified in RECRTS
IWK* 	integer work space of length at least LL (o/p from RECRTS)

      SUBROUTINE DOSCAL(E,A,B,LEN,DOS,WK,NW)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LEN),B(LEN),WK(NW,4)

Computes the density of states corresponding to a given continued fraction using the Christoffel function formalism of Haydock.

It also outputs the 2 kinds of real polynomial solutions to the symmetric recurrence, and the Christoffel function from which the result is calculated.
 
The two solutions to the symmetric recurrence:
	B(n+1)R(n+1) = (E-A(n))R(n) ö B(n)R(n-1), n=0,1,..,(LEN-1)
are initialised with
	P(0)=1, P(-1)=0 and Q(0)=0, B(0)Q(-1)=-1

The approximation to the DOS for i levels is then given by CFTL(i)/(pi*sqrt(inv.cur(i))) for i=1 to LEN

The output is ready for input to DYCFTL if the rate of decay of the Christoffel function is required.

ARGUMENTS : (* denotes overwritten argument)

      E	argument of the continued fraction
      A	diagonal coefficients of the cf (L-1 values)
      B	off-diagonal coefficients of the cf (L values)
      LEN	length of the cf
      DOS*	calculated value of the density of states
      WK*	work array.  On output contains:
 	Col 1	the values of the Christoffel function for each level
	Col 2	the values of the second function inv.cur for each level
	Col 3	the values of the first kind polynomials P(n)
	Col 4	the values of the second kind polynomials Q(n)
      NW	first DIMENSION of array WK in the calling segment






















      SUBROUTINE DYCFTL(LEN,LFIT,CFTL,EMAT,NE,EDIAG,DRATE,WK)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(LEN),B(LEN),CFTL(LEN),EMAT(NE,2),EDIAG(NE),WK(LFIT)

Estimates the decay rate of the Christoffel function after the latter has been evaluated by DOSCAL.  It is assumed to have a power-law behaviour, which is approximated by a linear least squares fit to the log of the function over the last LFIT values.  The last LFIT values of the tabulated function are used in the approximation.  If the argument of the cf was outside the support of the density of states, the decay rate being exponential, this is likely to give an unphysical large negative value.

The routine BS is used, and the routine PREPAS must be called first to factorise the rectangular matrix using HH.

ARGUMENTS : (* indicates an overwritten argument)

      LEN	length of continued fraction
      LFIT	number of levels used to estimated the decay rate
      CFTL	output from the routine DOSCAL, called with corresponding arguments. (The first col of the WK array.)  The last LFIT values are to be approximated.
      EMAT	unchanged output of PREPAS
      NE	first DIMENSION of array EMAT in the calling segment
      EDIAG	unchanged output of PREPAS
      DRATE*	estimated decay rate of the Christoffel function
      WK*	the O/P from BS.  The sum of squares of the last LFIT-2 values gives the sum of squares of deviations in the fit

      SUBROUTINE GRASSZ(Z,A,B,LEN,GRVAL,WK,NW)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LEN),B(LEN)
      COMPLEX*16 Z,GRVAL(2),WK(NW,2)

Evaluates the value of the Greenian corresponding to a given J-fraction (3-term recurrence relation) for a complex argument.  This is done by direction evaluation of the (symmetric) recurrence for the 2 polynomial solutions and applying the boundary condition that the Îtailâ is converged.

The two solutions to the symmetric recurrence:
	B(n+1)R(n+1) = (E-A(n))R(n) ö B(n)R(n-1), n=0,1,..,(LEN-1)
are initialised with
	P(0)=1/B(1), P(-1)=0 and Q(0)=0, Q(-1)=-1

ARGUMENTS : (* denotes overwritten argument)

      Z	complex argument
      A	diagonal elements of the recurrence
      B	off-diagonal elements of the recurrence
      LEN	length of the recurrence
      GRVAL*	complex estimated values of the Greenian on each sheet
      WK*	complex work-space containing the values of the polynomials on 
	output; the first kind in the first column.
      NW	first dimension of arrays WK and WKR


      SUBROUTINE GRCALZ(Z,A,B,LEN,GRVAL,WK,WKR,NW)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LEN),B(LEN),WKR(NW,2)
      COMPLEX*16 Z,GRVAL,WK(NW,3)

      Evaluates the value of the Greenian corresponding to a given J-fraction (3-term recurrence relation) on the Îphysicalâ sheet, for a complex argument.  This is done by direction evaluation of the (symmetric) recurrence for the 2 polynomial solutions, inverse summation for the sums of squares and direct summation for their pairwise product.
The two solutions to the symmetric recurrence:
	B(n+1)R(n+1) = (E-A(n))R(n) ö B(n)R(n-1), n=0,1,..,(LEN-1)
are initialised with
	P(0)=1/B(1), P(-1)=0 and Q(0)=0, Q(-1)=-1

ARGUMENTS : (* denotes overwritten argument)

      Z	complex argument
      A	diagonal elements of the recurrence
      B	off-diagonal elements of the recurrence
      LEN	length of the recurrence
      GRVAL*	complex estimated value of the Greenian
       WK*	complex work-space on output:
	Col 1	the values of the first kind polynomials P(n)
	Col 2	the values of the second kind polynomials Q(n)
	Col 3	the values of successive approximations to the Greenian
       WKR*	real work-space on output:
	Col 1	the values of the Christoffel function for each level for P
	Col 2	the values of the Christoffel function for each level for Q
       NW	first dimension of arrays WK and WKR


























      SUBROUTINE GRGENZ(Z,A,B,LEN,GRVAL,WK,WKR,NW)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(LEN),B(LEN),WKR(NW,3)
      COMPLEX*16 Z,GRVAL,WK(NW,2)

Evaluates the value of the Greenian corresponding to a given J-fraction (3-term recurrence relation) on the Îphysicalâ sheet, for a complex argument.  This is done by direction evaluation of the (symmetric) recurrence for the 2 polynomial solutions, direct double summation for their outer product, and inverse summation for the sums of squares.
The two solutions to the symmetric recurrence:
	B(n+1)R(n+1) = (E-A(n))R(n) ö B(n)R(n-1), n=0,1,..,(LEN-1)
are initialised with
	P(0)=1, P(-1)=0 and Q(0)=0, Q(-1)=-1/B(1)

ARGUMENTS : (* denotes overwritten argument)

      Z	complex argument
      A	diagonal elements of the recurrence
      B	off-diagonal elements of the recurrence
      LEN	length of the recurrence
      GRVAL*	complex estimated value of the Greenian
      WK*	complex work-space containing the values of the polynomials on 
	output; the first kind in the first column.
       WKR*	real work-space containing on output in the
	first column the Christoffel function of the first polynomials
	second column the inverse outer product of the two kinds
	third column the inner product of the two kinds
      NW	first dimension of arrays WK and WKR






      SUBROUTINE HH(A,MM,M,N,PT,EMACH)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(MM,N),PT(N)
 
Factorises the rectangular matrix A into an upper triangular matrix
and an orthogonal matrix using the Householder transformation without
interchanges.  The solution (in the least squares sense) to the
equation  Ax=b  or  (Atrans)Ax=(Atrans)b  may then be calculated
using BS
 
ARGUMENTS : (* denotes overwritten argument)
 
      A*	on input, the matrix to be factorised
      	on output, the upper triangle contains the transformed matrix,
      	and the lower triangle contains transformation information
      MM	first DIMENSION of array A in calling routine
      M	number of rows in matrix A
      N	 number of columns in matrix A (.le.M)
      PT*	the diagonal elements of the triangular matrix
      EMACH	machine accuracy

SUBROUTINE INFGEN(IC,N,WT,RNS,WTS,NPTS)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
              DIMENSION RNS(N),WTS(N)

Generates a quadrature for integrals over the infinite real line of functions multiplied by weight function .  This is based on the method of QFINT transformed to the infinite range using x = t/sqrt(1-t*t).  The number of points generated will be 1 less than a power of  2.

ARGUMENTS : (* indicates an overwritten argument)

      IC	code to indicate weight function required:
	4         gives exp(-x**2) (Gaussian weight)
	5         gives 1/(exp(pi*x/2)+exp(-pi*x/2))  (Meixner weight)
	6         gives exp(-x) for x>0, 0 for x<0 (Laguerre weight)
      N	maximum number of nodes required
      WT	total weight required
      RNS*	nodes
      WTS*	weights
      NPTS*	number of nodes generated  (= 2**n  - 1)



      SUBROUTINE INITAL(MMAX)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLKI/XV(511),WV(511)
 
Tabulates nodes and weights for quadrature using QFINT etc. in the order in which they are needed . The nodes used are actually XV(I)*(B-A)+A. The number of points generated is actually 1 less than integer power of 2.The weights used in the M point quadrature are { the first M values in WV divided by M+1} *(B-A) .  (M is 1 less than an integer power of two)
 
ARGUMENTS:
MMAX	Maximum number of function evaluations in any subsequent call to QFINT etc. <512 .




      FUNCTION NUMC(A,B2,ALAM,LM1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LM1),B2(LM1)
 
     Evaluates the number of eigenvalues of a symmetric tridiagonal matrix strictly greater than ALAM. The Sturm sequence property is used. No sub-diagonal element may be zero.[Wilkinson]
 
Arguments:
 
     NUMC takes the integer value required
     A    diagonal matrix elements i=1,LM1
     B2   squares of the subdiagonal matrix elements i=2,LM1
     ALAM point of evaluation
     LM1  dimension of matrix

      SUBROUTINE NUMD(A,B2,ALAM,LM1,DE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LM1),B2(LM1)
 
 
         Evaluates the inverse of the logarithmic derivative of the determinant (symmetric tridiagonal matrix - ALAM * identity) when no sub-diagonal matrix element is zero.[Wilkinson]
 
Arguments : ( * indicates an overwritten argument)
 
     A	diagonal matrix elements i=1,LM1
     B2	squares of sub-diagonal matrix elements i=2,LM1
     ALAM	argument of determinant above
     LM1	dimension of matrix
     DE*	det(matrix - ALAM)/det'(matrix - ALAM)

















      SUBROUTINE PREPAS(LEN,LFIT,EMACH,EMAT,NE,EDIAG)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION EMAT(NE,2),EDIAG(NE)

Prepare matrix for the least squares fitting to the log of a power law decay  of the Christoffel function by the routine DYCFTL.  The last LFIT values of the LEN given will be used in this.

This routine calls HH.

ARGUMENTS : (* indicates an overwritten argument)

      LEN	length of the recurrence provided
      LFIT	number of values to be used in the fitting
      EMACH	machine accuracy
      EMAT*	matrix ready for input to DYCFTL
      NE	first DIMENSION of array EMAT in the calling segment
      EDIAG*	vector ready for input to DYCFTL





      SUBROUTINE QFINT(A,B,F,FINT,ERR,E,MMAX,IP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLKI/XV(511),WV(511)
 
Evaluates the definite integral of a function F(X) between limits A and B . The method is a modified Gauss-Chebyshev (see C.M.M.Nex M.Sc. dissertation ) and does not call for function evaluations at the ends of the range .The first call of QFINT in a program must be preceded by a call to INITAL which sets up tables in the COMMON block/BLKI/ .
 
Copies of this routine are filed as QFINT1 and QFINT2 so that multiple  integrals may be evaluated by repeated application of the one-dimensional rule.
 
ARGUMENTS : (* indicates argument overwritten by the subroutine)
 
A	lower limit of range of integration
B	upper limit of range of integration
F	name of a REAL FUNCTION segment defining the integrand , with
	one REAL argument . It must be declared in an EXTERNAL statement
	at the head of the calling routine .
FINT*	the computed integral
ERR	 the required maximum absolute error
E*	the estimate of the error in the answer
MMAX	maximum number of function evaluations allowed <512 for current
	dimensions . If more are needed either use QQFINT or increase
	/BLKI/ dimensions in QFINT and INITAL to a suitable power of 2
IP	control integer set to the stream number to which  intermediate
	results  may be sent  ; set to 0 if not required.

























      SUBROUTINE QVINT(ENDS,N,F,FINT,ERR,E,MMAX,IP,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ENDS(N,2),WK(N)
      COMMON /BLKI/XV(511),WV(511)

Evaluates the definite integral of a function F(X) along a line in N-space . The method is a modified Gauss-Chebyshev (see C.M.M.Nex M.Sc. dissertation ) and does not call for function evaluations at the ends of the range.The first call of QVINT in a program must be preceded by a call to INITAL which sets up tables in the COMMON block/BLKI/ .
 
ARGUMENTS : (* indicates argument overwritten by the subroutine)
 
ENDS	the first column defines the lower limit of range of integration
	the second column defines the upper limit of range of integration
N	dimension of the space in which the geometry is defined
F	name of a REAL FUNCTION segment defining the integrand, with two arguments (Z,N)
	Z   a REAL vector defining the point  at which the function is to be evaluated
	N   DIMENSION of the first argument of F
	F must be declared in an EXTERNAL statement at the head of the calling routine
FINT*	the computed integral
ERR	the required maximum absolute error
E*	the estimate of the error in the answer
MMAX	maximum number of function evaluations allowed <512 for current dimensions.
	If more are needed increase /BLKI/ dimensions in QVINT and INITAL to a suitable power of 2
IP	control integer set to the stream number to which  intermediate results  may be sent; 
	set to 0 if no output required
WK*	work array used to communicate with F




























      SUBROUTINE QVINT1(CORN,N,F,FINT,ERR,E,MMAX,IP,AVAC,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CORN(N,3),WK(N,2)
      COMMON /BLKI/XV(511),WV(511)
      EXTERNAL F

Evaluates the definite integral of a function F(X) on a parallelogram in N-space . The method is a modified Gauss-Chebyshev (see C.M.M.Nex M.Sc. dissertation ) and does not call for function evaluations at the ends of the range.The first call of QVINT1 in a program must be preceded by a call to INITAL which sets up tables in the COMMON block/BLKI/ .
 
This routine calls FVINT1 and QVINT

ARGUMENTS : (* indicates argument overwritten by the subroutine)
 
CORN	the first column defines origin (bottom left corner) of the parallelogram
	the second column defines the lower right corner of the parallelogram (x)
	the third column defines the upper left corner of the parallelogram (y)
N	dimension of the space in which the geometry is defined
F	name of a REAL FUNCTION segment defining the integrand, with two arguments (Z,N)
	Z   a REAL vector defining the point at which the function is to be evaluated
	N   DIMENSION of the first argument of F
	F must be declared in an EXTERNAL statement at the head of the calling routine
FINT*	the computed integral
ERR	the required maximum absolute error
E*	the estimate of the error in the answer
MMAX	maximum number of function evaluations allowed per dimension <512 for current dimensions .
	If more are needed increase /BLKI/ dimensions in QVINT, QVINT1 and INITAL to a suitable power of 2
IP	control integer set to the stream number to which  intermediate results  may be sent; 
	set to 0 if no output required
AVAC*	average estimated error in the inner integral
WK*	work array used to communicate with FVINT1

      FUNCTION FVINT1(ARG,F,CORN,N,ERR,MMAX,ACC1,IC1)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARG(N,2),CORN(N,4)
      EXTERNAL F

Used by QVINT1 to set up the inner integral for a double integral















      SUBROUTINE QVINT2(CORN,N,F,FINT,ERR,E,MMAX,IP,AVAC,WK)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION CORN(N,4),WK(N,3),AVAC(2)
      COMMON /BLKI/XV(511),WV(511)
      EXTERNAL F


Evaluates the definite integral of a function F(X) in a parallelepiped in N-space. The method is a modified Gauss-Chebyshev (see C.M.M.Nex M.Sc. dissertation ) and does not call for function evaluations at the ends of the range. The first call of QVINT2 in a program must be preceded by a call to INITAL which sets up tables in the COMMON block/BLKI/.
 
This routine calls FVINT1, FVINT2, QVINT and QVINT1

ARGUMENTS : (* indicates argument overwritten by the subroutine)
 
CORN	the first column defines origin (bottom left corner) of the parallelepiped
	the second column defines the lower right corner of the parallelepiped (x)
	the third column defines the back left corner of the parallelepiped (y)
	the forth column defines the upper corner of the parallelepiped above the origin (z)
N	dimension of the space in which the geometry is defined
F	name of a REAL FUNCTION segment defining the integrand,
	with two arguments (Z,N)
	Z   a REAL vector defining the point  at which the function is to be evaluated
	N   DIMENSION of the first argument of F
	F must be declared in an EXTERNAL statement at the head of the calling routine
FINT*	the computed integral
ERR	the required maximum absolute error
E*	the estimate of the error in the answer
MMAX	maximum number of function evaluations allowed per dimension <512 for current dimensions.
	If more are needed increase /BLKI/ dimensions in QVINT, QVINT1, QVINT2 and INITAL to a suitable power of 2
IP	control integer set to the stream number to which  intermediate results  may be sent; 
	set to 0 if no output required
AVAC*	average estimated error in the two inner integrals:
	AVAC(1) corresponds to the innermost integral
	AVAC(2) corresponds to the middle integral
WK*	work array used to communicate with FVINT2


      FUNCTION FVINT2(ARG,F,CORN,N,ERR,MMAX,VACC,ICNT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARG(N,3),CORN(N,4),VACC(2)
      EXTERNAL F

Used by QVINT2 to set up the inner two integrals for the triple integral.









      SUBROUTINE RECQD(A,B2,LM1,X,WT,M,EPS,WK,MU)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LM1),B2(LM1),X(LM1),WT(LM1),WK(LM1),MU(LM1)
 
     Calculates the coefficients of the Gaussian quadrature corresponding to a given J-matrix or three-term recurrence relation. The Sturm sequence property is used to locate the nodes of the quadrature and the expression given in RECWTS to calculate the weights. [Cheney]
 
This routine calls RECWT, RECRTS, NUMC & NUMD
 
Arguments : (* indicates an overwritten argument)
 
     A	diagonal elements of the J-matrix i=1,LM1
     B2	squares of the off-diagonal elements of the J-matrix i=2,LM1
     	B2(1) gives the normalisation of the quadrature
     LM1	dimension of the J-matrix
     X*	nodes of Gaussian quadrature
     WT*	weights of Gaussian quadrature
     M*	number of quadrature nodes . If different from LM1 then
     	the routine has faulted.
     EPS	accuracy required in node calculation
     WK*	work array of length at least LM1
     MU*	work array of length at least LM1 (o/p from RECRTS)




      SUBROUTINE RECRTS(A,B2,LM1,EPS,XLIM,N,X,MULT,BI,NI)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LM1),B2(LM1),X(LM1),MULT(LM1),BI(LM1),NI(LM1)
 
 
     Computes some or all of the eigenvalues of a symmetric tridiagonal matrix with NO zero sub-diagonal elements (i.e.B2(i)>0 ).  The method used is bisection based on the Sturm sequence property followed by Newton's method for isolated roots.[Wilkinson]
 
This routine uses NUMC, NUMD
 
Arguments : (* indicates an overwritten argument)
 
     A	diagonal matrix elements i=1,LM1
     B2	squares of sub-diagonal matrix elements i=2,LM1
     LM1	dimension of matrix
     EPS	absolute accuracy required in eigenvalues
     XLIM	upper bound on eigenvalues to be found (if relevant)
     N*	if 0 on input : only eigenvalues less than XLIM are found
	on output : number of distinct eigenvalues found
     X*	eigenvalues in ascending order
     MULT	multiplicity of each eigenvalue
	if negative then the corresponding eigenvalue was found with less accuracy than EPS
     BI*	real work array of length at least LM1
     NI*	integer work array of length at least LM1

      SUBROUTINE RECSTD(A,B,IC,WT,N,EPS,RNS,WTS,WK,NWK)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION RNS(N),WTS(N),WK(NWK,4)

Generates a summation inner product corresponding to a given standard weight function, scaled and shifted as required.  This is made from the 3-term recurrence for IC=1 or 5, directly for IC=2 or 3, and via QFINT quadrature for IC=4,5,6.  In the latter cases it is advisable to use e.g. N=255,511 if the intention is to generate 20 levels of a continued fraction.  If the cf coefficients are generated they are in the monic form.

This routine calls RECQD, RECWT, RECRTS, NUMC, NUMD, INFGEN

ARGUMENTS:
    A	shift of standard distribution
    B	width renormalisation of standard distribution
    IC*	INPUT code for type of weight
		0    for delta function
		1    for constant (top hat)
		2    for inverse root(1-x*x)
		3    for root (1-x*x)
		4    for exp(-x*x)
		5    for Meixner weight
		6    for Laguerre weight (exp(-x) for x positive)
	    7   for Pollaczek weight
	    8    for Charlier weight
	OUTPUT -1 for failure, otherwise unchanged
		failure is only possible for non-Cheb weights
    WT	total weight or normalisation of the band
    N*	number of terms required in the summation
	overwritten with the number actually generated
    EPS	accuracy to be used in RECQD etc
    RNS*	nodes in summation
    WTS*	weights in summation
    WK*	work-space.
		On output the first two columns contain the cf cfts.
		when Meixner, Legendre, Pollaczek, Charlier (IC=1,5,6,7) used
    NWK	first DIMENSION of array WK in calling segment
	Greater than or equal to N

















      FUNCTION RECWT(E,A,B2,LL,EPS,N,P,NS)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(LL),B2(LL),P(2,3)

     Computes the value of the weight at the fixed point in a 1-fixed point Gaussian quadrature, given the corresponding 3-term recurrence relation:

          p(e,j)= (e-A(j)) * p(e,j-1) - B2(j) * p(e,j-2)

This routine may be called repeatedly with increasing number of levels such that it does not recompute earlier polynomial values.  If required the value of the last coefficient ,A(LL), may be computed, or it may be assumed that this has already been done and that value used in the calculation of the weight.
     The expression for the weight used is (with n=LL)
 
          { p(E,n-1)*q(E,n-1) - p(E,n)*q(E,n-2)} /
     { p(E,n-1)*p'(E,n) - p'(E,n-1)*p(E,n) + p(E,n)**2/B2(n) }
 
as this form is independent of the normalisation of the polynomials. p and q are the monic polynomials of the first and second kinds.
 
Arguments : (* indicates an overwritten argument)
 
RECWT     takes the value of the required weight
     E	required fixed point in quadrature. It may be a node of
	the LL-1 or LL quadrature if A(LL) is appropriately defined
     A*	diagonal elements of the recurrence. If N is changed
	from -1 input to 0 on output then A(LL) contains the adjusted
	value to achieve a Gaussian node at E, otherwise A is unchanged.
     B2	sub-diagonal elements of the recurrence
     LL	index of last B2 value to be used
     EPS	relative threshold value of the polynomial below which E will
	be accepted as a zero
     N*	code :
	 -1   A(LL) to be overwritten. N changed to 0 if successful,
	unchanged otherwise
	0   A(LL) given (not overwritten)
	1   A(LL) not computed explicitly (not overwritten)
     P*	final polynomial values used in calculation of weight
	to be used unchanged if routine is re-entered with NS=LL
	If n=LL-IABS(N)
               P(2,1)=p(E,n)       P(1,1)=p(E,n-1)
               P(2,2)=p'(E,n)      P(1,2)=p'(E,n-1)
               P(2,3)=q(E,n-1)     P(1,3)=q(E,n-2)
	q(E,m) is the polynomial of the second kind of degree m
     NS	point at which recurrence is initiated . This should be 1 initially , 
	but for a subsequent call, with E unchanged and larger LL, should
	be set to the current value of LL








      SUBROUTINE SIMBIS(FUN,A,B,EPS,NIT,X)
      IMPLICIT REAL*8(A-H,O-Z)

Uses simple bisection to locate a root of a function in a given interval.  If there is no sign difference between the ends of the interval, and error return is made

Arguments : (* indicates an overwritten argument)

FUN	name of an double precision FUNCTION segment with one double precision argument, defining the function whose zero is to be located.
	This must be declared EXTERNAL in the calling segment.
A	lower limit of range to be searched
B	upper limit of range to be searched
EPS	absolute accuracy required in the root position
NIT*	output code:
	0	no change of sign in the interval
	otherwise the number of iterations used
X*	value of the root







      SUBROUTINE VECCF(Z,AV,BETA,NV,NVE,LEN,TERM,VAL)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Z(NVE),AV(NV,LEN),BETA(LEN),TERM(NVE),VAL(NVE)

Evaluates a vector continued fraction as defined in Haydock, Nex & Wexler
The dimension of cf is one greater than dimension the diagonal vector coefficients, and conjugation of the vector is negation of the last (extended) component.

Arguments : (* indicates an overwritten argument)

Z	vector argument of the c.f.
A	diagonal vector parameters of the c.f.
BETA	scalar off-diagonal parameters of the c.f.
NV	dimension of the vector c.f. parameters, and 
	first DIMENSION of array AV
NVE	NV+1 : dimension of the argument of the c.f
	DIMENSION of arrays Z,TERM,VAL
LEN	length of the continued fraction
TERM	value of the vector terminator at the current point
VAL*	output value of the vector c.f.











     SUBROUTINE VECPOL(AV,BETA,NV,LEN,ARG,PV,PSCAL,CFTL)
     IMPLICIT REAL*8 (A-H,O-Z)
     DIMENSION AV(NV,LEN),BETA(LEN),ARG(NV),
    1          PV(NV,LEN),PSCAL(LEN),CFTL(LEN)

Evaluate in the Îreal subspaceâ the symmetric vector polynomial recurrence and 
accumulate the Christoffel function.

Arguments : (* indicates an overwritten argument)

AV	diagonal vector c.f. coefficients, first subscript giving the component index and the second the level index.
BETA	off-diagonal vector c.f.coefficients
NV	dimension of the vector c.f. and first DIMENSION of AV
LEN	length of the continued fraction
ARG	vector argument of the c.f.
PV*	vector 'polynomials' for each level
PSCAL*	scalar 'polynomials' for each level
CFTL*	Christoffel function for each level