c To get dgamma, "send dgamma from fnlib". c To get d1mach, mail netlib c send d1mach from core c subroutine gaussq(kind,n,alpha,beta,kpts,endpts,b,g,w) implicit none include 'max.inc' c c this set of routines computes the nodes t(j) and weights c w(j) for gaussian-type quadrature rules with pre-assigned c nodes. these are used when one wishes to approximate c c integral (from a to b) f(x) w(x) dx c c n c by sum w f(t ) c j=1 j j c c (note w(x) and w(j) have no connection with each other.) c here w(x) is one of six possible non-negative weight c functions (listed below), and f(x) is the c function to be integrated. gaussian quadrature is particularly c useful on infinite intervals (with appropriate weight c functions), since then other techniques often fail. c c associated with each weight function w(x) is a set of c orthogonal polynomials. the nodes t(j) are just the zeroes c of the proper n-th degree polynomial. c c input parameters (all real numbers are in double precision) c c kind an integer between 1 and 6 giving the type of c quadrature rule: c c kind = 1: legendre quadrature, w(x) = 1 on (-1, 1) c kind = 2: chebyshev quadrature of the first kind c w(x) = 1/sqrt(1 - x*x) on (-1, +1) c kind = 3: chebyshev quadrature of the second kind c w(x) = sqrt(1 - x*x) on (-1, 1) c kind = 4: hermite quadrature, w(x) = exp(-x*x) on c (-infinity, +infinity) c kind = 5: jacobi quadrature, w(x) = (1-x)**alpha * (1+x)** c beta on (-1, 1), alpha, beta .gt. -1. c note: kind=2 and 3 are a special case of this. c kind = 6: generalized laguerre quadrature, w(x) = exp(-x)* c x**alpha on (0, +infinity), alpha .gt. -1 c c n the number of points used for the quadrature rule c alpha real parameter used only for gauss-jacobi and gauss- c laguerre quadrature (otherwise use 0.d0). c beta real parameter used only for gauss-jacobi quadrature-- c (otherwise use 0.d0) c kpts (integer) normally 0, unless the left or right end- c point (or both) of the interval is required to be a c node (this is called gauss-radau or gauss-lobatto c quadrature). then kpts is the number of fixed c endpoints (1 or 2). c endpts real array of length 2. contains the values of c any fixed endpoints, if kpts = 1 or 2. c b real scratch array of length n c c output parameters (both double precision arrays of length n) c c t will contain the desired nodes. c w will contain the desired weights w(j). c c underflow may sometimes occur, but is harmless. c c references c 1. golub, g. h., and welsch, j. h., "calculation of gaussian c quadrature rules," mathematics of computation 23 (april, c 1969), pp. 221-230. c 2. golub, g. h., "some modified matrix eigenvalue problems," c siam review 15 (april, 1973), pp. 318-334 (section 7). c 3. stroud and secrest, gaussian quadrature formulas, prentice- c hall, englewood cliffs, n.j., 1966. c c original version 20 jan 1975 from stanford c modified 21 dec 1983 by eric grosse c imtql2 => gausq2 c hex constant => d1mach (from core library) c compute pi using datan c removed accuracy claims, description of method c added single precision version c c inputs integer kind integer n double precision alpha double precision beta integer kpts double precision endpts(1:2) double precision b(1:Nqmx) c outputs double precision g(1:Nqmx) double precision w(1:Nqmx) c temp double precision t(1:Nqmx) double precision w1(1:Nqmx) double precision a(1:Nmax,1:Nc_mx) double precision sum_w integer i,ierr double precision muzero double precision t1,gam,solve,dsqrt c label integer strlen character*(Nchar_mx) label label='subroutine gaussq' c Debug write(*,*) 'kind=',kind c Debug c call class(kind,n,alpha,beta,b,t,muzero) c c the matrix of coefficients is assumed to be symmetric. c the array t contains the diagonal elements, the array c b the off-diagonal elements. c make appropriate changes in the lower right 2 by 2 c submatrix. c if (kpts.eq.0) then goto 100 endif if (kpts.eq.2) then goto 50 endif c c if kpts=1, only t(n) must be changed c t(n)=solve(endpts(1),n,t,b)*b(n-1)**2+endpts(1) goto 100 c c if kpts=2, t(n) and b(n-1) must be recomputed c 50 gam=solve(endpts(1),n,t,b) t1=((endpts(1)-endpts(2))/(solve(endpts(2),n,t,b)-gam)) b(n-1)=dsqrt(t1) t(n)=endpts(1)+gam*t1 c c note that the indices of the elements of b run from 1 to n-1 c and thus the value of b(n) is arbitrary. c now compute the eigenvalues of the symmetric tridiagonal c matrix, which has been modified as necessary. c the method used is a ql-type method with origin shifting c 100 w(1)=1.0D+0 do i=2,n w(i)=0.0D0 enddo c call gausq2(n,t,b,w,ierr) do i=1,n w(i)=muzero*w(i)*w(i) enddo c do i=1,n a(i,1)=t(i) a(i,2)=w(i) enddo call ssort(a,n,2) do i=1,n t(i)=a(i,1) w1(i)=a(i,2) enddo c modification of abscissas t(i) and weights w1(i) in order to c produce abscissas g(i) and weights w(i) sum_w=0.0D+0 do i=1,n if (kind.ne.0) then g(i)=(t(i)+1.0D+0)/2.0D+0 w(i)=w1(i)/2.0D+0 endif sum_w=sum_w+w(i) enddo if (dabs(sum_w-1.0D+0).gt.1.0D-5) then write(*,*) 'warning: sum of weights=',sum_w c stop endif return end c c c double precision function solve(shift,n,a,b) c c this procedure performs elimination to solve for the c n-th component of the solution delta to the equation c c (jn - shift*identity) * delta = en, c c where en is the vector of all zeroes except for 1 in c the n-th position. c c the matrix jn is symmetric tridiagonal, with diagonal c elements a(i), off-diagonal elements b(i). this equation c must be solved to obtain the appropriate changes in the lower c 2 by 2 submatrix of coefficients for orthogonal polynomials. c c implicit none include 'max.inc' integer n,nm1,i double precision shift,a(1:Nqmx),b(1:Nqmx),alpha c alpha=a(1)-shift nm1=n-1 do i=2,nm1 alpha=a(i)-shift-b(i-1)**2/alpha enddo solve=1.0D+0/alpha return end c c c subroutine class(kind,n,alpha,beta,b,a,muzero) implicit none c c this procedure supplies the coefficients a(j), b(j) of the c recurrence relation c c b p (x) = (x - a ) p (x) - b p (x) c j j j j-1 j-1 j-2 c c for the various classical (normalized) orthogonal polynomials, c and the zero-th moment c c muzero = integral w(x) dx c c of the given polynomial's weight function w(x). since the c polynomials are orthonormalized, the tridiagonal matrix is c guaranteed to be symmetric. c c the input parameter alpha is used only for laguerre and c jacobi polynomials, and the parameter beta is used only for c jacobi polynomials. the laguerre and jacobi polynomials c require the gamma function. c include 'max.inc' double precision a(1:Nqmx),b(1:Nqmx),muzero,alpha,beta double precision abi,a2b2,dgamma,pi,dsqrt,ab integer kind,n,nm1,i c pi=4.0D+0*datan(1.0D+0) nm1=n-1 go to (10, 20, 30, 40, 50, 60), kind c c kind = 1: legendre polynomials p(x) c on (-1, +1), w(x) = 1. c 10 muzero = 2.0d0 do 11 i = 1, nm1 a(i) = 0.0d0 abi = i 11 b(i) = abi/dsqrt(4*abi*abi - 1.0d0) a(n) = 0.0d0 return c c kind = 2: chebyshev polynomials of the first kind t(x) c on (-1, +1), w(x) = 1 / sqrt(1 - x*x) c 20 muzero = pi do 21 i = 1, nm1 a(i) = 0.0d0 21 b(i) = 0.5d0 b(1) = dsqrt(0.5d0) a(n) = 0.0d0 return c c kind = 3: chebyshev polynomials of the second kind u(x) c on (-1, +1), w(x) = sqrt(1 - x*x) c 30 muzero = pi/2.0d0 do 31 i = 1, nm1 a(i) = 0.0d0 31 b(i) = 0.5d0 a(n) = 0.0d0 return c c kind = 4: hermite polynomials h(x) on (-infinity, c +infinity), w(x) = exp(-x**2) c 40 muzero = dsqrt(pi) do 41 i = 1, nm1 a(i) = 0.0d0 41 b(i) = dsqrt(i/2.0d0) a(n) = 0.0d0 return c c kind = 5: jacobi polynomials p(alpha, beta)(x) on c (-1, +1), w(x) = (1-x)**alpha + (1+x)**beta, alpha and c beta greater than -1 c 50 ab = alpha + beta abi = 2.0d0 + ab muzero = 2.0d0 ** (ab + 1.0d0) * dgamma(alpha + 1.0d0) * dgamma( x beta + 1.0d0) / dgamma(abi) a(1) = (beta - alpha)/abi b(1) = dsqrt(4.0d0*(1.0d0 + alpha)*(1.0d0 + beta)/((abi + 1.0d0)* 1 abi*abi)) a2b2 = beta*beta - alpha*alpha do 51 i = 2, nm1 abi = 2.0d0*i + ab a(i) = a2b2/((abi - 2.0d0)*abi) 51 b(i) = dsqrt (4.0d0*i*(i + alpha)*(i + beta)*(i + ab)/ 1 ((abi*abi - 1)*abi*abi)) abi = 2.0d0*n + ab a(n) = a2b2/((abi - 2.0d0)*abi) return c c kind = 6: laguerre polynomials l(alpha)(x) on c (0, +infinity), w(x) = exp(-x) * x**alpha, alpha greater c than -1. c 60 muzero = dgamma(alpha + 1.0d0) do 61 i = 1, nm1 a(i) = 2.0d0*i - 1.0d0 + alpha 61 b(i) = dsqrt(i*(i + alpha)) a(n) = 2.0d0*n - 1 + alpha return end c c subroutine gausq2(n,d,e,z,ierr) implicit none c c this subroutine is a translation of an algol procedure, c num. math. 12, 377-383(1968) by martin and wilkinson, c as modified in num. math. 15, 450(1970) by dubrulle. c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). c this is a modified version of the 'eispack' routine imtql2. c c this subroutine finds the eigenvalues and first components of the c eigenvectors of a symmetric tridiagonal matrix by the implicit ql c method. c c on input: c c n is the order of the matrix; c c d contains the diagonal elements of the input matrix; c c e contains the subdiagonal elements of the input matrix c in its first n-1 positions. e(n) is arbitrary; c c z contains the first row of the identity matrix. c c on output: c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1, 2, ..., ierr-1; c c e has been destroyed; c c z contains the first components of the orthonormal eigenvectors c of the symmetric tridiagonal matrix. if an error exit is c made, z contains the eigenvectors associated with the stored c eigenvalues; c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c ------------------------------------------------------------------ c include 'max.inc' integer i, j, k, l, m, n, ii, mml, ierr real*8 d(1:Nqmx),e(1:Nqmx),z(1:Nqmx), b, c, f, g, p, r, s, machep real*8 dsqrt, dabs, dsign, d1mach c machep=d1mach(4) c ierr = 0 if (n .eq. 1) go to 1001 c e(n) = 0.0d0 do 240 l = 1, n j = 0 c :::::::::: look for small sub-diagonal element :::::::::: 105 do 110 m = l, n if (m .eq. n) go to 120 if (dabs(e(m)) .le. machep * (dabs(d(m)) + dabs(d(m+1)))) x go to 120 110 continue c 120 p = d(l) if (m .eq. l) go to 240 if (j .eq. 30) go to 1000 j = j + 1 c :::::::::: form shift :::::::::: g = (d(l+1) - p) / (2.0d0 * e(l)) r = dsqrt(g*g+1.0d0) g = d(m) - p + e(l) / (g + dsign(r, g)) s = 1.0d0 c = 1.0d0 p = 0.0d0 mml = m - l c c :::::::::: for i=m-1 step -1 until l do -- :::::::::: do 200 ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) if (dabs(f) .lt. dabs(g)) go to 150 c = g / f r = dsqrt(c*c+1.0d0) e(i+1) = f * r s = 1.0d0 / r c = c * s go to 160 150 s = f / g r = dsqrt(s*s+1.0d0) e(i+1) = g * r c = 1.0d0 / r s = s * c 160 g = d(i+1) - p r = (d(i) - g) * s + 2.0d0 * c * b p = s * r d(i+1) = g + p g = c * r - b c :::::::::: form first component of vector :::::::::: f = z(i+1) z(i+1) = s * z(i) + c * f 200 z(i) = c * z(i) - s * f c d(l) = d(l) - p e(l) = g e(m) = 0.0d0 go to 105 240 continue c c :::::::::: order eigenvalues and eigenvectors :::::::::: do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p p = z(i) z(i) = z(k) z(k) = p 300 continue c go to 1001 c :::::::::: set error -- no convergence to an c eigenvalue after 30 iterations :::::::::: 1000 ierr = l 1001 return c :::::::::: last card of gausq2 :::::::::: end c c c double precision function dgamma(x) double precision x dgamma = 1.0d0 return end *DECK D1MACH DOUBLE PRECISION FUNCTION D1MACH (I) C***BEGIN PROLOGUE D1MACH C***PURPOSE Return floating point machine dependent constants. C***LIBRARY SLATEC C***CATEGORY R1 C***TYPE DOUBLE PRECISION (R1MACH-S, D1MACH-D) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Fox, P. A., (Bell Labs) C Hall, A. D., (Bell Labs) C Schryer, N. L., (Bell Labs) C***DESCRIPTION C C D1MACH can be used to obtain machine-dependent parameters for the C local machine environment. It is a function subprogram with one C (input) argument, and can be referenced as follows: C C D = D1MACH(I) C C where I=1,...,5. The (output) value of D above is determined by C the (input) value of I. The results for various values of I are C discussed below. C C D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. C D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. C D1MACH( 3) = B**(-T), the smallest relative spacing. C D1MACH( 4) = B**(1-T), the largest relative spacing. C D1MACH( 5) = LOG10(B) C C Assume double precision numbers are represented in the T-digit, C base-B form C C sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) C C where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and C EMIN .LE. E .LE. EMAX. C C The values of B, T, EMIN and EMAX are provided in I1MACH as C follows: C I1MACH(10) = B, the base. C I1MACH(14) = T, the number of base-B digits. C I1MACH(15) = EMIN, the smallest exponent E. C I1MACH(16) = EMAX, the largest exponent E. C C To alter this function for a particular environment, the desired C set of DATA statements should be activated by removing the C from C column 1. Also, the values of D1MACH(1) - D1MACH(4) should be C checked for consistency with the local operating system. C C***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for C a portable library, ACM Transactions on Mathematical C Software 4, 2 (June 1978), pp. 177-188. C***ROUTINES CALLED XERMSG C***REVISION HISTORY (YYMMDD) C 750101 DATE WRITTEN C 890213 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ) C 900618 Added DEC RISC constants. (WRB) C 900723 Added IBM RS 6000 constants. (WRB) C 900911 Added SUN 386i constants. (WRB) C 910710 Added HP 730 constants. (SMR) C 911114 Added Convex IEEE constants. (WRB) C 920121 Added SUN -r8 compiler option constants. (WRB) C 920229 Added Touchstone Delta i860 constants. (WRB) C 920501 Reformatted the REFERENCES section. (WRB) C 920625 Added CONVEX -p8 and -pd8 compiler option constants. C (BKS, WRB) C 930201 Added DEC Alpha and SGI constants. (RWC and WRB) C***END PROLOGUE D1MACH C INTEGER SMALL(4) INTEGER LARGE(4) INTEGER RIGHT(4) INTEGER DIVER(4) INTEGER LOG10(4) C DOUBLE PRECISION DMACH(5) SAVE DMACH C EQUIVALENCE (DMACH(1),SMALL(1)) EQUIVALENCE (DMACH(2),LARGE(1)) EQUIVALENCE (DMACH(3),RIGHT(1)) EQUIVALENCE (DMACH(4),DIVER(1)) EQUIVALENCE (DMACH(5),LOG10(1)) C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING THE 68020/68881 COMPILER OPTION C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE AMIGA C ABSOFT FORTRAN COMPILER USING SOFTWARE FLOATING POINT C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FDFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE APOLLO C C DATA SMALL(1), SMALL(2) / 16#00100000, 16#00000000 / C DATA LARGE(1), LARGE(2) / 16#7FFFFFFF, 16#FFFFFFFF / C DATA RIGHT(1), RIGHT(2) / 16#3CA00000, 16#00000000 / C DATA DIVER(1), DIVER(2) / 16#3CB00000, 16#00000000 / C DATA LOG10(1), LOG10(2) / 16#3FD34413, 16#509F79FF / C C MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM C C DATA SMALL(1) / ZC00800000 / C DATA SMALL(2) / Z000000000 / C DATA LARGE(1) / ZDFFFFFFFF / C DATA LARGE(2) / ZFFFFFFFFF / C DATA RIGHT(1) / ZCC5800000 / C DATA RIGHT(2) / Z000000000 / C DATA DIVER(1) / ZCC6800000 / C DATA DIVER(2) / Z000000000 / C DATA LOG10(1) / ZD00E730E7 / C DATA LOG10(2) / ZC77800DC0 / C C MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O0000000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O0007777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS C C DATA SMALL(1) / O1771000000000000 / C DATA SMALL(2) / O7770000000000000 / C DATA LARGE(1) / O0777777777777777 / C DATA LARGE(2) / O7777777777777777 / C DATA RIGHT(1) / O1461000000000000 / C DATA RIGHT(2) / O0000000000000000 / C DATA DIVER(1) / O1451000000000000 / C DATA DIVER(2) / O0000000000000000 / C DATA LOG10(1) / O1157163034761674 / C DATA LOG10(2) / O0006677466732724 / C C MACHINE CONSTANTS FOR THE CDC 170/180 SERIES USING NOS/VE C C DATA SMALL(1) / Z"3001800000000000" / C DATA SMALL(2) / Z"3001000000000000" / C DATA LARGE(1) / Z"4FFEFFFFFFFFFFFE" / C DATA LARGE(2) / Z"4FFE000000000000" / C DATA RIGHT(1) / Z"3FD2800000000000" / C DATA RIGHT(2) / Z"3FD2000000000000" / C DATA DIVER(1) / Z"3FD3800000000000" / C DATA DIVER(2) / Z"3FD3000000000000" / C DATA LOG10(1) / Z"3FFF9A209A84FBCF" / C DATA LOG10(2) / Z"3FFFF7988F8959AC" / C C MACHINE CONSTANTS FOR THE CDC 6000/7000 SERIES C C DATA SMALL(1) / 00564000000000000000B / C DATA SMALL(2) / 00000000000000000000B / C DATA LARGE(1) / 37757777777777777777B / C DATA LARGE(2) / 37157777777777777777B / C DATA RIGHT(1) / 15624000000000000000B / C DATA RIGHT(2) / 00000000000000000000B / C DATA DIVER(1) / 15634000000000000000B / C DATA DIVER(2) / 00000000000000000000B / C DATA LOG10(1) / 17164642023241175717B / C DATA LOG10(2) / 16367571421742254654B / C C MACHINE CONSTANTS FOR THE CELERITY C1260 C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fn OR -pd8 COMPILER OPTION C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CC0000000000000' / C DATA DMACH(4) / Z'3CD0000000000000' / C DATA DMACH(5) / Z'3FF34413509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -fi COMPILER OPTION C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE CONVEX C USING THE -p8 COMPILER OPTION C C DATA DMACH(1) / Z'00010000000000000000000000000000' / C DATA DMACH(2) / Z'7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3F900000000000000000000000000000' / C DATA DMACH(4) / Z'3F910000000000000000000000000000' / C DATA DMACH(5) / Z'3FFF34413509F79FEF311F12B35816F9' / C C MACHINE CONSTANTS FOR THE CRAY C C DATA SMALL(1) / 201354000000000000000B / C DATA SMALL(2) / 000000000000000000000B / C DATA LARGE(1) / 577767777777777777777B / C DATA LARGE(2) / 000007777777777777774B / C DATA RIGHT(1) / 376434000000000000000B / C DATA RIGHT(2) / 000000000000000000000B / C DATA DIVER(1) / 376444000000000000000B / C DATA DIVER(2) / 000000000000000000000B / C DATA LOG10(1) / 377774642023241175717B / C DATA LOG10(2) / 000007571421742254654B / C C MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 C NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING CARD - C STATIC DMACH(5) C C DATA SMALL / 20K, 3*0 / C DATA LARGE / 77777K, 3*177777K / C DATA RIGHT / 31420K, 3*0 / C DATA DIVER / 32020K, 3*0 / C DATA LOG10 / 40423K, 42023K, 50237K, 74776K / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING G_FLOAT C C DATA DMACH(1) / '0000000000000010'X / C DATA DMACH(2) / 'FFFFFFFFFFFF7FFF'X / C DATA DMACH(3) / '0000000000003CC0'X / C DATA DMACH(4) / '0000000000003CD0'X / C DATA DMACH(5) / '79FF509F44133FF3'X / C C MACHINE CONSTANTS FOR THE DEC ALPHA C USING IEEE_FORMAT C C DATA DMACH(1) / '0010000000000000'X / C DATA DMACH(2) / '7FEFFFFFFFFFFFFF'X / C DATA DMACH(3) / '3CA0000000000000'X / C DATA DMACH(4) / '3CB0000000000000'X / C DATA DMACH(5) / '3FD34413509F79FF'X / C C MACHINE CONSTANTS FOR THE DEC RISC C C DATA SMALL(1), SMALL(2) / Z'00000000', Z'00100000'/ C DATA LARGE(1), LARGE(2) / Z'FFFFFFFF', Z'7FEFFFFF'/ C DATA RIGHT(1), RIGHT(2) / Z'00000000', Z'3CA00000'/ C DATA DIVER(1), DIVER(2) / Z'00000000', Z'3CB00000'/ C DATA LOG10(1), LOG10(2) / Z'509F79FF', Z'3FD34413'/ C C MACHINE CONSTANTS FOR THE DEC VAX C USING D_FLOATING C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA LOG10(1), LOG10(2) / 546979738, -805796613 / C C DATA SMALL(1), SMALL(2) / Z00000080, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00002480, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00002500, Z00000000 / C DATA LOG10(1), LOG10(2) / Z209A3F9A, ZCFF884FB / C C MACHINE CONSTANTS FOR THE DEC VAX C USING G_FLOATING C (EXPRESSED IN INTEGER AND HEXADECIMAL) C THE HEX FORMAT BELOW MAY NOT BE SUITABLE FOR UNIX SYSTEMS C THE INTEGER FORMAT SHOULD BE OK FOR UNIX SYSTEMS C C DATA SMALL(1), SMALL(2) / 16, 0 / C DATA LARGE(1), LARGE(2) / -32769, -1 / C DATA RIGHT(1), RIGHT(2) / 15552, 0 / C DATA DIVER(1), DIVER(2) / 15568, 0 / C DATA LOG10(1), LOG10(2) / 1142112243, 2046775455 / C C DATA SMALL(1), SMALL(2) / Z00000010, Z00000000 / C DATA LARGE(1), LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z00003CC0, Z00000000 / C DATA DIVER(1), DIVER(2) / Z00003CD0, Z00000000 / C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / C DATA RIGHT(1), RIGHT(2) / '3CB00000'X,'00000000'X / C DATA DIVER(1), DIVER(2) / '3CC00000'X,'00000000'X / C DATA LOG10(1), LOG10(2) / '3FD34413'X,'509F79FF'X / C C MACHINE CONSTANTS FOR THE HARRIS 220 C C DATA SMALL(1), SMALL(2) / '20000000, '00000201 / C DATA LARGE(1), LARGE(2) / '37777777, '37777577 / C DATA RIGHT(1), RIGHT(2) / '20000000, '00000333 / C DATA DIVER(1), DIVER(2) / '20000000, '00000334 / C DATA LOG10(1), LOG10(2) / '23210115, '10237777 / C C MACHINE CONSTANTS FOR THE HONEYWELL 600/6000 SERIES C C DATA SMALL(1), SMALL(2) / O402400000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O376777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O604400000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O606400000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O776464202324, O117571775714 / C C MACHINE CONSTANTS FOR THE HP 730 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE HP 2100 C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / C DATA RIGHT(1), RIGHT(2), RIGHT(3) / 40000B, 0, 265B / C DATA DIVER(1), DIVER(2), DIVER(3) / 40000B, 0, 276B / C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / C DATA LARGE(1), LARGE(2) / 77777B, 177777B / C DATA LARGE(3), LARGE(4) / 177777B, 177776B / C DATA RIGHT(1), RIGHT(2) / 40000B, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 225B / C DATA DIVER(1), DIVER(2) / 40000B, 0 / C DATA DIVER(3), DIVER(4) / 0, 227B / C DATA LOG10(1), LOG10(2) / 46420B, 46502B / C DATA LOG10(3), LOG10(4) / 76747B, 176377B / C C MACHINE CONSTANTS FOR THE HP 9000 C C DATA SMALL(1), SMALL(2) / 00040000000B, 00000000000B / C DATA LARGE(1), LARGE(2) / 17737777777B, 37777777777B / C DATA RIGHT(1), RIGHT(2) / 07454000000B, 00000000000B / C DATA DIVER(1), DIVER(2) / 07460000000B, 00000000000B / C DATA LOG10(1), LOG10(2) / 07764642023B, 12047674777B / C C MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, C THE XEROX SIGMA 5/7/9, THE SEL SYSTEMS 85/86, AND C THE PERKIN ELMER (INTERDATA) 7/32. C C DATA SMALL(1), SMALL(2) / Z00100000, Z00000000 / C DATA LARGE(1), LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / C DATA RIGHT(1), RIGHT(2) / Z33100000, Z00000000 / C DATA DIVER(1), DIVER(2) / Z34100000, Z00000000 / C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE IBM PC C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. C C DATA SMALL(1) / 2.23D-308 / C DATA LARGE(1) / 1.79D+308 / C DATA RIGHT(1) / 1.11D-16 / C DATA DIVER(1) / 2.22D-16 / C DATA LOG10(1) / 0.301029995663981195D0 / C C MACHINE CONSTANTS FOR THE IBM RS 6000 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE INTEL i860 C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR) C C DATA SMALL(1), SMALL(2) / "033400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "344777777777 / C DATA RIGHT(1), RIGHT(2) / "113400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "114400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "144117571776 / C C MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR) C C DATA SMALL(1), SMALL(2) / "000400000000, "000000000000 / C DATA LARGE(1), LARGE(2) / "377777777777, "377777777777 / C DATA RIGHT(1), RIGHT(2) / "103400000000, "000000000000 / C DATA DIVER(1), DIVER(2) / "104400000000, "000000000000 / C DATA LOG10(1), LOG10(2) / "177464202324, "476747767461 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 8388608, 0 / C DATA LARGE(1), LARGE(2) / 2147483647, -1 / C DATA RIGHT(1), RIGHT(2) / 612368384, 0 / C DATA DIVER(1), DIVER(2) / 620756992, 0 / C DATA LOG10(1), LOG10(2) / 1067065498, -2063872008 / C C DATA SMALL(1), SMALL(2) / O00040000000, O00000000000 / C DATA LARGE(1), LARGE(2) / O17777777777, O37777777777 / C DATA RIGHT(1), RIGHT(2) / O04440000000, O00000000000 / C DATA DIVER(1), DIVER(2) / O04500000000, O00000000000 / C DATA LOG10(1), LOG10(2) / O07746420232, O20476747770 / C C MACHINE CONSTANTS FOR PDP-11 FORTRAN SUPPORTING C 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). C C DATA SMALL(1), SMALL(2) / 128, 0 / C DATA SMALL(3), SMALL(4) / 0, 0 / C DATA LARGE(1), LARGE(2) / 32767, -1 / C DATA LARGE(3), LARGE(4) / -1, -1 / C DATA RIGHT(1), RIGHT(2) / 9344, 0 / C DATA RIGHT(3), RIGHT(4) / 0, 0 / C DATA DIVER(1), DIVER(2) / 9472, 0 / C DATA DIVER(3), DIVER(4) / 0, 0 / C DATA LOG10(1), LOG10(2) / 16282, 8346 / C DATA LOG10(3), LOG10(4) / -31493, -12296 / C C DATA SMALL(1), SMALL(2) / O000200, O000000 / C DATA SMALL(3), SMALL(4) / O000000, O000000 / C DATA LARGE(1), LARGE(2) / O077777, O177777 / C DATA LARGE(3), LARGE(4) / O177777, O177777 / C DATA RIGHT(1), RIGHT(2) / O022200, O000000 / C DATA RIGHT(3), RIGHT(4) / O000000, O000000 / C DATA DIVER(1), DIVER(2) / O022400, O000000 / C DATA DIVER(3), DIVER(4) / O000000, O000000 / C DATA LOG10(1), LOG10(2) / O037632, O020232 / C DATA LOG10(3), LOG10(4) / O102373, O147770 / C C MACHINE CONSTANTS FOR THE SILICON GRAPHICS C C DATA SMALL(1), SMALL(2) / Z'00100000', Z'00000000' / C DATA LARGE(1), LARGE(2) / Z'7FEFFFFF', Z'FFFFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'3CA00000', Z'00000000' / C DATA DIVER(1), DIVER(2) / Z'3CB00000', Z'00000000' / C DATA LOG10(1), LOG10(2) / Z'3FD34413', Z'509F79FF' / C C MACHINE CONSTANTS FOR THE SUN C C DATA DMACH(1) / Z'0010000000000000' / C DATA DMACH(2) / Z'7FEFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3CA0000000000000' / C DATA DMACH(4) / Z'3CB0000000000000' / C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE SUN C USING THE -r8 COMPILER OPTION C C DATA DMACH(1) / Z'00010000000000000000000000000000' / C DATA DMACH(2) / Z'7FFEFFFFFFFFFFFFFFFFFFFFFFFFFFFF' / C DATA DMACH(3) / Z'3F8E0000000000000000000000000000' / C DATA DMACH(4) / Z'3F8F0000000000000000000000000000' / C DATA DMACH(5) / Z'3FFD34413509F79FEF311F12B35816F9' / C C MACHINE CONSTANTS FOR THE SUN 386i C C DATA SMALL(1), SMALL(2) / Z'FFFFFFFD', Z'000FFFFF' / C DATA LARGE(1), LARGE(2) / Z'FFFFFFB0', Z'7FEFFFFF' / C DATA RIGHT(1), RIGHT(2) / Z'000000B0', Z'3CA00000' / C DATA DIVER(1), DIVER(2) / Z'FFFFFFCB', Z'3CAFFFFF' C DATA LOG10(1), LOG10(2) / Z'509F79E9', Z'3FD34413' / C C MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES FTN COMPILER C C DATA SMALL(1), SMALL(2) / O000040000000, O000000000000 / C DATA LARGE(1), LARGE(2) / O377777777777, O777777777777 / C DATA RIGHT(1), RIGHT(2) / O170540000000, O000000000000 / C DATA DIVER(1), DIVER(2) / O170640000000, O000000000000 / C DATA LOG10(1), LOG10(2) / O177746420232, O411757177572 / C C***FIRST EXECUTABLE STATEMENT D1MACH IF (I .LT. 1 .OR. I .GT. 5) then write(*,*) 'error in d1mach' c CALL XERMSG ('SLATEC', 'D1MACH', c + 'I OUT OF BOUNDS', 1, 2) endif C D1MACH = DMACH(I) RETURN C END