!THIS FILE CONTAINS SUBROUTINES TO PERFORM THE FOLLOWING OPERATIONS: ! ! SUBROUTINE PRODUCES !---------------------------------- ------------------------------------ ! ! 1 GETQUAD2(QUAD,DG,N,VIEW_FLAG, VECTORS OF QUADRATURE ROOTS AND ! VIEW_ANGLE,MU,W) WEIGHTS FOR GAUSS, DOUBLE GAUSS, ! OR LOBATTO QUADRATURE WITH COSINE ! OF USER ANGLE & 0 WEIGHT APPENDED ! TO THE RESPECTIVE VECTORS ! (USES NETLIB SUBS) ! 2 INVERT(N,A,AINV,INV_IO) THE INVERSE OF A MATRIX ! (USES LINPACK SUBS) ! 3 MATDIAG(N,VEC,DIAG) AN NxN DIAGONAL MATRIX ! 4 MATIDENT(N,IDENTITY) AN NxN IDENTITY MATRIX ! 5 MM_DG(N,DVEC,G,PROD) THE PRODUCT OF A DIAGONAL MATRIX AND ! A GENERAL MATRIX ! 6 MM_D1GD2(N,D1VEC,G,D2VEC,PROD) THE PRODUCT OF A GENERAL MATRIX AND ! TWO DIAGONAL MATRICES ! 7 MM_GD(N,G,DVEC,PROD) THE PRODUCT OF A GENERAL MATRIX AND ! A DIAGONAL MATRIX ! 8 MM_IG1G2(N,M,G1,G2,PROD,INV_IO) THE PRODUCT OF THE INVERSE OF A ! GENERAL MATRIX G1 AND A GENERAL ! MATRIX G2 ! (USES LINPACK SUBS) ! 9 MV_DV(N,DVEC,V,PROD) THE PRODUCT OF A DIAGONAL MATRIX AND ! A VECTOR !10 PLEG(ANGLEO,GMU,M,NEXP,NZEN, THE RENORMALIZED ASSOCIATED LEGENDRE ! YPLEG) POLYNOMIALS !11 SYSSOL(N,M,A,B,X,SYS_IO) THE SOLUTION OF A LINEAR SYSTEM OF ! EQUATIONS VIA GAUSSIAN ELIMINATION ! (USES LINPACK SUBS) ! !************************************************************************ !************************************************************************ SUBROUTINE GETQUAD2(QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE,MU,W) !INPUT : QUAD,DG,N,VIEW_FLAG,VIEW_ANGLE !OUTPUT: MU,W !THIS PROGRAM OBTAINS THE MUs AND WEIGHTS FOR GAUSSIAN OR LABATTO !QUADRATURE NEEDED FOR RADIANT'S RT COMPUTATIONS. !PROGRAMMER: MATT CHRISTI !DATE: 5/1/03 !DATA DICTIONARY******************************************************* ! ! DG = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL ! BE USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED ! (0=NORMAL GAUSS,1=DOUBLE GAUSS) ! MU = VECTOR HOLDING QUADRATURE ROOTS (COSINES OF QUADRATURE ! ANGLES) ! N = NUMBER OF UPWARD (OR EQUIVALENTLY DOWNWARD) STREAMS ! NS = TOTAL NUMBER OF STREAMS ! QUAD = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED IN ! COMPUTING THE VALUES OF MU AND W (0=GAUSS,1=LABATTO) ! VIEW_ANGLE = THE USER-DEFINED VIEW ANGLE (in degrees from zenith) ! VIEW_FLAG = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE ! (0=OFF, 1=ON) ! W = VECTOR HOLDING QUADRATURE WEIGHTS ! !********************************************************************** !INTRINSIC SUBPROGRAMS USED BY GETQUAD2******************************** ! NONE !********************************************************************** !EXTERNAL SUBPROGRAMS USED BY GETQUAD2********************************* ! QDATA !********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER :: QUAD,DG,N,VIEW_FLAG DOUBLE PRECISION :: VIEW_ANGLE !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N) :: MU,W !INTERNAL VARIABLES INTEGER :: J,NS,KPTS,Q_IO,QUAD_IO,VEC_SIZE DOUBLE PRECISION :: PI PARAMETER (PI=3.1415926535897932D0) DOUBLE PRECISION :: VIEW_MU DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: & MUTEMP,WTEMP !PROGRAM START QUAD_IO = 0 IF(QUAD_IO == 1) THEN PRINT* PRINT*,'ENTERING GETQUAD2' END IF NS=2*N IF (QUAD == 0) THEN KPTS = QUAD ELSE IF (QUAD == 1) THEN KPTS = QUAD + 1 END IF Q_IO = 0 IF (VIEW_FLAG == 1) THEN !USING USER-DEFINED ANGLE & PRE-DEFINED QUADRATURE ANGLES IF (DG == 1) THEN !USING "DOUBLE-GAUSS" VEC_SIZE = N-1 ELSE !USING NORMAL GAUSS VEC_SIZE = NS-2 END IF ELSE !USING PRE-DEFINED QUADRATURE ANGLES ONLY IF (DG == 1) THEN !USING "DOUBLE-GAUSS" VEC_SIZE = N ELSE !USING NORMAL GAUSS VEC_SIZE = NS END IF END IF ALLOCATE( MUTEMP(VEC_SIZE),WTEMP(VEC_SIZE) ) CALL QDATA(VEC_SIZE,KPTS,MUTEMP,WTEMP,Q_IO) MUTEMP = -1.0D0*MUTEMP IF (DG == 1) THEN !USING "DOUBLE-GAUSS" IF(QUAD_IO == 1) THEN PRINT*,'BEFORE DG:' DO J=1,VEC_SIZE WRITE(*,5) J,MUTEMP(J),WTEMP(J) 5 FORMAT(2X,I2,2(3X,E15.9E2)) END DO END IF !SWITCHING TO "DOUBLE GAUSS" DO J=1,VEC_SIZE MUTEMP(J) = 0.5D0*MUTEMP(J) + 0.5D0 WTEMP(J) = 0.5D0*WTEMP(J) END DO END IF IF(QUAD_IO == 1) THEN DO J=1,VEC_SIZE WRITE(*,5) J,MUTEMP(J),WTEMP(J) END DO END IF IF (VIEW_FLAG == 1) THEN VIEW_MU = DCOS(VIEW_ANGLE/180.0D0*PI) DO J=1,N IF (J == N) THEN MU(J) = VIEW_MU W(J) = 0.0D0 ELSE MU(J) = MUTEMP(J) W(J) = WTEMP(J) END IF END DO ELSE DO J=1,N MU(J) = MUTEMP(J) W(J) = WTEMP(J) END DO END IF IF(QUAD_IO == 1) THEN PRINT* DO J=1,N WRITE(*,5) J,MU(J),W(J) END DO END IF DEALLOCATE (MUTEMP,WTEMP) IF(QUAD_IO == 1) THEN PRINT* PRINT*,'LEAVING GETQUAD2' END IF RETURN END !********************************************************************** !********************************************************************** SUBROUTINE QDATA(N,KPTS,ROOT,WEIGHT,Q_IO) !THIS SUBROUTINE PRESETS MOST THE PARAMETERS FOR THE GAUSSQ SUBROUTINE ! TO RETURN THE QUADRATURE ROOTS AND WEIGHTS FOR THE LEGENDRE ! POLYNOMIALS ! WRITTEN BY M. CHRISTI ! DATE: 08/28/01 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! TO USE: ! ! (1) SET N FOR TOTAL NUMBER OF STREAMS YOU ARE USING ! (2) SET KPTS AS BELOW FOR DESIRED QUADRATURE SCHEME: ! *FOR PURE GAUSSIAN QUADRATURE, SET KPTS = 0 ! *FOR GAUSS-LOBATTO QUADRATURE, SET KPTS = 2 ! (3) SET Q_IO TO 0 OR 1 (QDATA I/O FLAG) ! *FOR NORMAL OPERATION, SET TO 0 ! *TO DISPLAY OUTPUT OF GAUSSQ FOR TESTING, SET TO 1 ! (4) PROVIDE TWO ONE-DIMENSIONAL ARRAYS OF SIZE N IN YOUR ROUTINE ! TO HOLD ROOTS AND WEIGHTS ! ! FOR REFERENCE, THE SUBROUTINE GAUSSQ AND OTHER SUBROUTINES UPON WHICH ! IT DEPENDS CAN BE OBTAINED FROM THE NETLIB REPOSITORY AT THE ! FOLLOWING URL: www.netlib.org/cgi-bin/search.pl ! AND TYPING "gaussq" IN THE SEARCH BOX !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! IMPLICIT NONE INTEGER I,N,Q_IO,KIND,KPTS DOUBLE PRECISION ALPHA,BETA,ENDPTS(2),B(N),ROOT(N),WEIGHT(N) IF (Q_IO == 1) THEN PRINT* PRINT*,'ENTERING QDATA' END IF !INPUT TO SUBROUTINE !TYPE OF QUADRATURE RULE KIND = 1 !FOR GAUSS-JACOBI AND GAUSS-LAGUERRE QUADRATURE ONLY ALPHA = 0.0D0 !FOR GAUSS-JACOBI QUADRATURE ONLY BETA = 0.0D0 !NUMBER OF FIXED ENDPOINTS ! KPTS = 2 !SPECIFIES FIXED ENDPOINTS ENDPTS(1) = 1.0D0 ENDPTS(2) = -1.0D0 !SCRATCH VECTOR B OF LENGTH N B = 0.0D0 !ROOT WILL CONTAIN VECTOR OF N GAUSSIAN ROOTS (I.E. THE NODES OR ! ABSCISSAE) ON RETURN ROOT = B !WEIGHT WILL CONTAIN VECTOR OF N GAUSSIAN WEIGHTS ON RETURN WEIGHT = B !COMPUTE GAUSSIAN ROOTS AND WEIGHTS CALL GAUSSQ(KIND,N,ALPHA,BETA,KPTS,ENDPTS,B,ROOT,WEIGHT) !DISPLAY OUTPUT IF DESIRED IF (Q_IO == 1) THEN IF (KPTS == 0) THEN WRITE(*,*) WRITE(*,*) 'USING PURE GAUSSIAN QUADRATURE' ELSE IF (KPTS == 2) THEN WRITE(*,*) WRITE(*,*) 'USING GAUSS-LOBATTO QUADRATURE' END IF WRITE(*,*) WRITE(*,5) 'FOR A VALUE OF N =',N WRITE(*,*) WRITE(*,*) ' THE ROOTS LOOK LIKE: THE WEIGHTS LOOK LIKE: ' DO I=1,N WRITE(*,10) ROOT(I), WEIGHT(I) END DO 5 FORMAT(A19,1X,I2) 10 FORMAT(1X,E22.15,3X,E22.15) PRINT* PRINT*,'LEAVING QDATA' END IF RETURN END !*********************************************************************** !*********************************************************************** subroutine gaussq(kind, n, alpha, beta, kpts, endpts, b, t, w) ! ! this set of routines computes the nodes t(j) and weights ! w(j) for gaussian-type quadrature rules with pre-assigned ! nodes. these are used when one wishes to approximate ! ! integral (from a to b) f(x) w(x) dx ! ! n ! by sum w f(t ) ! j=1 j j ! ! (note w(x) and w(j) have no connection with each other.) ! here w(x) is one of six possible non-negative weight ! functions (listed below), and f(x) is the ! function to be integrated. gaussian quadrature is particularly ! useful on infinite intervals (with appropriate weight ! functions), since then other techniques often fail. ! ! associated with each weight function w(x) is a set of ! orthogonal polynomials. the nodes t(j) are just the zeroes ! of the proper n-th degree polynomial. ! ! input parameters (all real numbers are in double precision) ! ! kind an integer between 1 and 6 giving the type of ! quadrature rule: ! ! kind = 1: legendre quadrature, w(x) = 1 on (-1, 1) ! kind = 2: chebyshev quadrature of the first kind ! w(x) = 1/sqrt(1 - x*x) on (-1, +1) ! kind = 3: chebyshev quadrature of the second kind ! w(x) = sqrt(1 - x*x) on (-1, 1) ! kind = 4: hermite quadrature, w(x) = exp(-x*x) on ! (-infinity, +infinity) ! kind = 5: jacobi quadrature, w(x) = (1-x)**alpha * (1+x)** ! beta on (-1, 1), alpha, beta .gt. -1. ! note: kind=2 and 3 are a special case of this. ! kind = 6: generalized laguerre quadrature, w(x) = exp(-x)* ! x**alpha on (0, +infinity), alpha .gt. -1 ! ! n the number of points used for the quadrature rule ! alpha real parameter used only for gauss-jacobi and gauss- ! laguerre quadrature (otherwise use 0.d0). ! beta real parameter used only for gauss-jacobi quadrature-- ! (otherwise use 0.d0) ! kpts (integer) normally 0, unless the left or right end- ! point (or both) of the interval is required to be a ! node (this is called gauss-radau or gauss-lobatto ! quadrature). then kpts is the number of fixed ! endpoints (1 or 2). ! endpts real array of length 2. contains the values of ! any fixed endpoints, if kpts = 1 or 2. ! b real scratch array of length n ! ! output parameters (both double precision arrays of length n) ! ! t will contain the desired nodes. ! w will contain the desired weights w(j). ! ! underflow may sometimes occur, but is harmless. ! ! references ! 1. golub, g. h., and welsch, j. h., "calculation of gaussian ! quadrature rules," mathematics of computation 23 (april, ! 1969), pp. 221-230. ! 2. golub, g. h., "some modified matrix eigenvalue problems," ! siam review 15 (april, 1973), pp. 318-334 (section 7). ! 3. stroud and secrest, gaussian quadrature formulas, prentice- ! hall, englewood cliffs, n.j., 1966. ! ! original version 20 jan 1975 from stanford ! modified 21 dec 1983 by eric grosse ! imtql2 => gausq2 ! hex constant => d1mach (from core library) ! compute pi using datan ! removed accuracy claims, description of method ! added single precision version ! double precision b(n), t(n), w(n), endpts(2), muzero, t1, & gam, solve, dsqrt, alpha, beta ! call class (kind, n, alpha, beta, b, t, muzero) ! ! the matrix of coefficients is assumed to be symmetric. ! the array t contains the diagonal elements, the array ! b the off-diagonal elements. ! make appropriate changes in the lower right 2 by 2 ! submatrix. ! if (kpts.eq.0) go to 100 if (kpts.eq.2) go to 50 ! ! if kpts=1, only t(n) must be changed ! t(n) = solve(endpts(1), n, t, b)*b(n-1)**2 + endpts(1) go to 100 ! ! if kpts=2, t(n) and b(n-1) must be recomputed ! 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 ! ! note that the indices of the elements of b run from 1 to n-1 ! and thus the value of b(n) is arbitrary. ! now compute the eigenvalues of the symmetric tridiagonal ! matrix, which has been modified as necessary. ! the method used is a ql-type method with origin shifting ! 100 w(1) = 1.0d0 do 105 i = 2, n 105 w(i) = 0.0d0 ! call gausq2 (n, t, b, w, ierr) do 110 i = 1, n 110 w(i) = muzero * w(i) * w(i) ! return end !********************************************************************** !********************************************************************** double precision function solve(shift, n, a, b) ! ! this procedure performs elimination to solve for the ! n-th component of the solution delta to the equation ! ! (jn - shift*identity) * delta = en, ! ! where en is the vector of all zeroes except for 1 in ! the n-th position. ! ! the matrix jn is symmetric tridiagonal, with diagonal ! elements a(i), off-diagonal elements b(i). this equation ! must be solved to obtain the appropriate changes in the lower ! 2 by 2 submatrix of coefficients for orthogonal polynomials. ! ! double precision shift, a(n), b(n), alpha ! alpha = a(1) - shift nm1 = n - 1 do 10 i = 2, nm1 10 alpha = a(i) - shift - b(i-1)**2/alpha solve = 1.0d0/alpha return end !********************************************************************** !********************************************************************** subroutine class(kind, n, alpha, beta, b, a, muzero) ! ! this procedure supplies the coefficients a(j), b(j) of the ! recurrence relation ! ! b p (x) = (x - a ) p (x) - b p (x) ! j j j j-1 j-1 j-2 ! ! for the various classical (normalized) orthogonal polynomials, ! and the zero-th moment ! ! muzero = integral w(x) dx ! ! of the given polynomial's weight function w(x). since the ! polynomials are orthonormalized, the tridiagonal matrix is ! guaranteed to be symmetric. ! ! the input parameter alpha is used only for laguerre and ! jacobi polynomials, and the parameter beta is used only for ! jacobi polynomials. the laguerre and jacobi polynomials ! require the gamma function. ! double precision a(n), b(n), muzero, alpha, beta double precision abi, a2b2, dgamma, pi, dsqrt, ab ! pi = 4.0d0 * datan(1.0d0) nm1 = n - 1 go to (10, 20, 30, 40, 50, 60), kind ! ! kind = 1: legendre polynomials p(x) ! on (-1, +1), w(x) = 1. ! 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 ! ! kind = 2: chebyshev polynomials of the first kind t(x) ! on (-1, +1), w(x) = 1 / sqrt(1 - x*x) ! 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 ! ! kind = 3: chebyshev polynomials of the second kind u(x) ! on (-1, +1), w(x) = sqrt(1 - x*x) ! 30 muzero = pi/2.0d0 do 31 i = 1, nm1 a(i) = 0.0d0 31 b(i) = 0.5d0 a(n) = 0.0d0 return ! ! kind = 4: hermite polynomials h(x) on (-infinity, ! +infinity), w(x) = exp(-x**2) ! 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 ! ! kind = 5: jacobi polynomials p(alpha, beta)(x) on ! (-1, +1), w(x) = (1-x)**alpha + (1+x)**beta, alpha and ! beta greater than -1 ! 50 ab = alpha + beta abi = 2.0d0 + ab muzero = 2.0d0 ** (ab + 1.0d0) * dgamma(alpha + 1.0d0) * & dgamma(beta + 1.0d0) / dgamma(abi) a(1) = (beta - alpha)/abi b(1) = dsqrt(4.0d0*(1.0d0 + alpha)*(1.0d0 + beta)/ & ((abi + 1.0d0)* 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)/ & ((abi*abi - 1)*abi*abi)) abi = 2.0d0*n + ab a(n) = a2b2/((abi - 2.0d0)*abi) return ! ! kind = 6: laguerre polynomials l(alpha)(x) on ! (0, +infinity), w(x) = exp(-x) * x**alpha, alpha greater ! than -1. ! 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 !********************************************************************** !********************************************************************** subroutine gausq2(n, d, e, z, ierr) ! ! this subroutine is a translation of an algol procedure, ! num. math. 12, 377-383(1968) by martin and wilkinson, ! as modified in num. math. 15, 450(1970) by dubrulle. ! handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). ! this is a modified version of the 'eispack' routine imtql2. ! ! this subroutine finds the eigenvalues and first components of the ! eigenvectors of a symmetric tridiagonal matrix by the implicit ql ! method. ! ! on input: ! ! n is the order of the matrix; ! ! d contains the diagonal elements of the input matrix; ! ! e contains the subdiagonal elements of the input matrix ! in its first n-1 positions. e(n) is arbitrary; ! ! z contains the first row of the identity matrix. ! ! on output: ! ! d contains the eigenvalues in ascending order. if an ! error exit is made, the eigenvalues are correct but ! unordered for indices 1, 2, ..., ierr-1; ! ! e has been destroyed; ! ! z contains the first components of the orthonormal eigenvectors ! of the symmetric tridiagonal matrix. if an error exit is ! made, z contains the eigenvectors associated with the stored ! eigenvalues; ! ! ierr is set to ! zero for normal return, ! j if the j-th eigenvalue has not been ! determined after 30 iterations. ! ! ------------------------------------------------------------------ ! integer i, j, k, l, m, n, ii, mml, ierr real*8 d(n), e(n), z(n), b, c, f, g, p, r, s, machep real*8 dsqrt, dabs, dsign, d1mach ! machep=d1mach(4) ! ierr = 0 if (n .eq. 1) go to 1001 ! e(n) = 0.0d0 do 240 l = 1, n j = 0 ! :::::::::: 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)))) & go to 120 110 continue ! 120 p = d(l) if (m .eq. l) go to 240 if (j .eq. 30) go to 1000 j = j + 1 ! :::::::::: 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 ! ! :::::::::: 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 ! :::::::::: 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 ! d(l) = d(l) - p e(l) = g e(m) = 0.0d0 go to 105 240 continue ! ! :::::::::: order eigenvalues and eigenvectors :::::::::: do 300 ii = 2, n i = ii - 1 k = i p = d(i) ! do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue ! 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 ! go to 1001 ! :::::::::: set error -- no convergence to an ! eigenvalue after 30 iterations :::::::::: 1000 ierr = l 1001 return ! :::::::::: last card of gausq2 :::::::::: end !********************************************************************** !********************************************************************** DOUBLE PRECISION FUNCTION D1MACH(I) INTEGER I ! ! DOUBLE-PRECISION MACHINE CONSTANTS ! D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. ! D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. ! D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. ! D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. ! D1MACH( 5) = LOG10(B) ! INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) INTEGER SC, CRAY1(38), J COMMON /D9MACH/ CRAY1 SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC DOUBLE PRECISION DMACH(5) 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)) ! THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES. ! R1MACH CAN HANDLE AUTO-DOUBLE COMPILING, BUT THIS VERSION OF ! D1MACH DOES NOT, BECAUSE WE DO NOT HAVE QUAD CONSTANTS FOR ! MANY MACHINES YET. ! TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 ! ON THE NEXT LINE DATA SC/0/ ! AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. ! CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY ! mail netlib@research.bell-labs.com ! send old1mach from blas ! PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. ! ! MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. ! DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / ! DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / ! DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / ! DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / ! DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ ! ! MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING ! 32-BIT INTEGERS. ! DATA SMALL(1),SMALL(2) / 8388608, 0 / ! DATA LARGE(1),LARGE(2) / 2147483647, -1 / ! DATA RIGHT(1),RIGHT(2) / 612368384, 0 / ! DATA DIVER(1),DIVER(2) / 620756992, 0 / ! DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ ! ! MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. ! DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / ! DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / ! DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / ! DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / ! DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ ! ! ON FIRST CALL, IF NO DATA UNCOMMENTED, TEST MACHINE TYPES. IF (SC .NE. 987) THEN DMACH(1) = 1.D13 IF ( SMALL(1) .EQ. 1117925532 & .AND. SMALL(2) .EQ. -448790528) THEN ! *** IEEE BIG ENDIAN *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2146435071 LARGE(2) = -1 RIGHT(1) = 1017118720 RIGHT(2) = 0 DIVER(1) = 1018167296 DIVER(2) = 0 LOG10(1) = 1070810131 LOG10(2) = 1352628735 ELSE IF ( SMALL(2) .EQ. 1117925532 & .AND. SMALL(1) .EQ. -448790528) THEN ! *** IEEE LITTLE ENDIAN *** SMALL(2) = 1048576 SMALL(1) = 0 LARGE(2) = 2146435071 LARGE(1) = -1 RIGHT(2) = 1017118720 RIGHT(1) = 0 DIVER(2) = 1018167296 DIVER(1) = 0 LOG10(2) = 1070810131 LOG10(1) = 1352628735 ELSE IF ( SMALL(1) .EQ. -2065213935 & .AND. SMALL(2) .EQ. 10752) THEN ! *** VAX WITH D_FLOATING *** SMALL(1) = 128 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 9344 RIGHT(2) = 0 DIVER(1) = 9472 DIVER(2) = 0 LOG10(1) = 546979738 LOG10(2) = -805796613 ELSE IF ( SMALL(1) .EQ. 1267827943 & .AND. SMALL(2) .EQ. 704643072) THEN ! *** IBM MAINFRAME *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 856686592 RIGHT(2) = 0 DIVER(1) = 873463808 DIVER(2) = 0 LOG10(1) = 1091781651 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 1120022684 & .AND. SMALL(2) .EQ. -448790528) THEN ! *** CONVEX C-1 *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 1019215872 RIGHT(2) = 0 DIVER(1) = 1020264448 DIVER(2) = 0 LOG10(1) = 1072907283 LOG10(2) = 1352628735 ELSE IF ( SMALL(1) .EQ. 815547074 & .AND. SMALL(2) .EQ. 58688) THEN ! *** VAX G-FLOATING *** SMALL(1) = 16 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 15552 RIGHT(2) = 0 DIVER(1) = 15568 DIVER(2) = 0 LOG10(1) = 1142112243 LOG10(2) = 2046775455 ELSE DMACH(2) = 1.D27 + 1 DMACH(3) = 1.D27 LARGE(2) = LARGE(2) - RIGHT(2) IF (LARGE(2) .EQ. 64 .AND. SMALL(2) .EQ. 0) THEN CRAY1(1) = 67291416 DO 10 J = 1, 20 CRAY1(J+1) = CRAY1(J) + CRAY1(J) 10 CONTINUE CRAY1(22) = CRAY1(21) + 321322 DO 20 J = 22, 37 CRAY1(J+1) = CRAY1(J) + CRAY1(J) 20 CONTINUE IF (CRAY1(38) .EQ. SMALL(1)) THEN ! *** CRAY *** CALL I1MCRY(SMALL(1), J, 8285, 8388608, 0) SMALL(2) = 0 CALL I1MCRY(LARGE(1), J, 24574, 16777215, 16777215) CALL I1MCRY(LARGE(2), J, 0, 16777215, 16777214) CALL I1MCRY(RIGHT(1), J, 16291, 8388608, 0) RIGHT(2) = 0 CALL I1MCRY(DIVER(1), J, 16292, 8388608, 0) DIVER(2) = 0 CALL I1MCRY(LOG10(1), J, 16383, 10100890, 8715215) CALL I1MCRY(LOG10(2), J, 0, 16226447, 9001388) ELSE WRITE(*,9000) STOP 779 END IF ELSE WRITE(*,9000) STOP 779 END IF END IF SC = 987 END IF ! SANITY CHECK IF (DMACH(4) .GE. 1.0D0) STOP 778 IF (I .LT. 1 .OR. I .GT. 5) THEN WRITE(*,*) 'D1MACH(I): I =',I,' is out of bounds.' STOP END IF D1MACH = DMACH(I) RETURN 9000 FORMAT(/' Adjust D1MACH by uncommenting data statements'/ & ' appropriate for your machine.') ! /* Standard C source for D1MACH -- remove the * in column 1 */ !#include !#include !#include !double d1mach_(long *i) !{ ! switch(*i){ ! case 1: return DBL_MIN; ! case 2: return DBL_MAX; ! case 3: return DBL_EPSILON/FLT_RADIX; ! case 4: return DBL_EPSILON; ! case 5: return log10((double)FLT_RADIX); ! } ! fprintf(stderr, "invalid argument: d1mach(%ld)\n", *i); ! exit(1); return 0; /* some compilers demand return values */ !} END !********************************************************************** !********************************************************************** SUBROUTINE I1MCRY(A, A1, B, C, D) !*** SPECIAL COMPUTATION FOR OLD CRAY MACHINES **** INTEGER A, A1, B, C, D A1 = 16777216*B + C A = 16777216*A1 + D END INTEGER FUNCTION I1MACH(I) INTEGER I ! ! I1MACH( 1) = THE STANDARD INPUT UNIT. ! I1MACH( 2) = THE STANDARD OUTPUT UNIT. ! I1MACH( 3) = THE STANDARD PUNCH UNIT. ! I1MACH( 4) = THE STANDARD ERROR MESSAGE UNIT. ! I1MACH( 5) = THE NUMBER OF BITS PER INTEGER STORAGE UNIT. ! I1MACH( 6) = THE NUMBER OF CHARACTERS PER CHARACTER STORAGE UNIT. ! INTEGERS HAVE FORM SIGN ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) ) ! I1MACH( 7) = A, THE BASE. ! I1MACH( 8) = S, THE NUMBER OF BASE-A DIGITS. ! I1MACH( 9) = A**S - 1, THE LARGEST MAGNITUDE. ! FLOATS HAVE FORM SIGN (B**E)*( (X(1)/B) + ... + (X(T)/B**T) ) ! WHERE EMIN .LE. E .LE. EMAX. ! I1MACH(10) = B, THE BASE. ! SINGLE-PRECISION ! I1MACH(11) = T, THE NUMBER OF BASE-B DIGITS. ! I1MACH(12) = EMIN, THE SMALLEST EXPONENT E. ! I1MACH(13) = EMAX, THE LARGEST EXPONENT E. ! DOUBLE-PRECISION ! I1MACH(14) = T, THE NUMBER OF BASE-B DIGITS. ! I1MACH(15) = EMIN, THE SMALLEST EXPONENT E. ! I1MACH(16) = EMAX, THE LARGEST EXPONENT E. ! INTEGER IMACH(16), OUTPUT, SC, SMALL(2) SAVE IMACH, SC REAL RMACH EQUIVALENCE (IMACH(4),OUTPUT), (RMACH,SMALL(1)) INTEGER I3, J, K, T3E(3) DATA T3E(1) / 9777664 / DATA T3E(2) / 5323660 / DATA T3E(3) / 46980 / ! THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES, ! INCLUDING AUTO-DOUBLE COMPILERS. ! TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 ! ON THE NEXT LINE DATA SC/0/ ! AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. ! CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY ! mail netlib@research.bell-labs.com ! send old1mach from blas ! PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. ! ! MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. ! ! DATA IMACH( 1) / 5 / ! DATA IMACH( 2) / 6 / ! DATA IMACH( 3) / 43 / ! DATA IMACH( 4) / 6 / ! DATA IMACH( 5) / 36 / ! DATA IMACH( 6) / 4 / ! DATA IMACH( 7) / 2 / ! DATA IMACH( 8) / 35 / ! DATA IMACH( 9) / O377777777777 / ! DATA IMACH(10) / 2 / ! DATA IMACH(11) / 27 / ! DATA IMACH(12) / -127 / ! DATA IMACH(13) / 127 / ! DATA IMACH(14) / 63 / ! DATA IMACH(15) / -127 / ! DATA IMACH(16) / 127 /, SC/987/ ! ! MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING ! 32-BIT INTEGER ARITHMETIC. ! ! DATA IMACH( 1) / 5 / ! DATA IMACH( 2) / 6 / ! DATA IMACH( 3) / 7 / ! DATA IMACH( 4) / 6 / ! DATA IMACH( 5) / 32 / ! DATA IMACH( 6) / 4 / ! DATA IMACH( 7) / 2 / ! DATA IMACH( 8) / 31 / ! DATA IMACH( 9) / 2147483647 / ! DATA IMACH(10) / 2 / ! DATA IMACH(11) / 24 / ! DATA IMACH(12) / -127 / ! DATA IMACH(13) / 127 / ! DATA IMACH(14) / 56 / ! DATA IMACH(15) / -127 / ! DATA IMACH(16) / 127 /, SC/987/ ! ! MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. ! ! NOTE THAT THE PUNCH UNIT, I1MACH(3), HAS BEEN SET TO 7 ! WHICH IS APPROPRIATE FOR THE UNIVAC-FOR SYSTEM. ! IF YOU HAVE THE UNIVAC-FTN SYSTEM, SET IT TO 1. ! ! DATA IMACH( 1) / 5 / ! DATA IMACH( 2) / 6 / ! DATA IMACH( 3) / 7 / ! DATA IMACH( 4) / 6 / ! DATA IMACH( 5) / 36 / ! DATA IMACH( 6) / 6 / ! DATA IMACH( 7) / 2 / ! DATA IMACH( 8) / 35 / ! DATA IMACH( 9) / O377777777777 / ! DATA IMACH(10) / 2 / ! DATA IMACH(11) / 27 / ! DATA IMACH(12) / -128 / ! DATA IMACH(13) / 127 / ! DATA IMACH(14) / 60 / ! DATA IMACH(15) /-1024 / ! DATA IMACH(16) / 1023 /, SC/987/ ! IF (SC .NE. 987) THEN ! *** CHECK FOR AUTODOUBLE *** SMALL(2) = 0 RMACH = 1E13 IF (SMALL(2) .NE. 0) THEN ! *** AUTODOUBLED *** IF ( (SMALL(1) .EQ. 1117925532 & .AND. SMALL(2) .EQ. -448790528) & .OR. (SMALL(2) .EQ. 1117925532 & .AND. SMALL(1) .EQ. -448790528)) THEN ! *** IEEE *** IMACH(10) = 2 IMACH(14) = 53 IMACH(15) = -1021 IMACH(16) = 1024 ELSE IF ( SMALL(1) .EQ. -2065213935 & .AND. SMALL(2) .EQ. 10752) THEN ! *** VAX WITH D_FLOATING *** IMACH(10) = 2 IMACH(14) = 56 IMACH(15) = -127 IMACH(16) = 127 ELSE IF ( SMALL(1) .EQ. 1267827943 & .AND. SMALL(2) .EQ. 704643072) THEN ! *** IBM MAINFRAME *** IMACH(10) = 16 IMACH(14) = 14 IMACH(15) = -64 IMACH(16) = 63 ELSE WRITE(*,9010) STOP 777 END IF IMACH(11) = IMACH(14) IMACH(12) = IMACH(15) IMACH(13) = IMACH(16) ELSE RMACH = 1234567. IF (SMALL(1) .EQ. 1234613304) THEN ! *** IEEE *** IMACH(10) = 2 IMACH(11) = 24 IMACH(12) = -125 IMACH(13) = 128 IMACH(14) = 53 IMACH(15) = -1021 IMACH(16) = 1024 SC = 987 ELSE IF (SMALL(1) .EQ. -1271379306) THEN ! *** VAX *** IMACH(10) = 2 IMACH(11) = 24 IMACH(12) = -127 IMACH(13) = 127 IMACH(14) = 56 IMACH(15) = -127 IMACH(16) = 127 SC = 987 ELSE IF (SMALL(1) .EQ. 1175639687) THEN ! *** IBM MAINFRAME *** IMACH(10) = 16 IMACH(11) = 6 IMACH(12) = -64 IMACH(13) = 63 IMACH(14) = 14 IMACH(15) = -64 IMACH(16) = 63 SC = 987 ELSE IF (SMALL(1) .EQ. 1251390520) THEN ! *** CONVEX C-1 *** IMACH(10) = 2 IMACH(11) = 24 IMACH(12) = -128 IMACH(13) = 127 IMACH(14) = 53 IMACH(15) = -1024 IMACH(16) = 1023 ELSE DO 10 I3 = 1, 3 J = SMALL(1) / 10000000 K = SMALL(1) - 10000000*J IF (K .NE. T3E(I3)) GO TO 20 SMALL(1) = J 10 CONTINUE ! *** CRAY T3E *** IMACH( 1) = 5 IMACH( 2) = 6 IMACH( 3) = 0 IMACH( 4) = 0 IMACH( 5) = 64 IMACH( 6) = 8 IMACH( 7) = 2 IMACH( 8) = 63 CALL I1MCR1(IMACH(9), K, 32767, 16777215, 16777215) IMACH(10) = 2 IMACH(11) = 53 IMACH(12) = -1021 IMACH(13) = 1024 IMACH(14) = 53 IMACH(15) = -1021 IMACH(16) = 1024 GO TO 35 20 CALL I1MCR1(J, K, 16405, 9876536, 0) IF (SMALL(1) .NE. J) THEN WRITE(*,9020) STOP 777 END IF ! *** CRAY 1, XMP, 2, AND 3 *** IMACH(1) = 5 IMACH(2) = 6 IMACH(3) = 102 IMACH(4) = 6 IMACH(5) = 46 IMACH(6) = 8 IMACH(7) = 2 IMACH(8) = 45 CALL I1MCR1(IMACH(9), K, 0, 4194303, 16777215) IMACH(10) = 2 IMACH(11) = 47 IMACH(12) = -8188 IMACH(13) = 8189 IMACH(14) = 94 IMACH(15) = -8141 IMACH(16) = 8189 GO TO 35 END IF END IF IMACH( 1) = 5 IMACH( 2) = 6 IMACH( 3) = 7 IMACH( 4) = 6 IMACH( 5) = 32 IMACH( 6) = 4 IMACH( 7) = 2 IMACH( 8) = 31 IMACH( 9) = 2147483647 35 SC = 987 END IF 9010 FORMAT(/' Adjust autodoubled I1MACH by uncommenting data'/ & ' statements appropriate for your machine and setting'/ & ' IMACH(I) = IMACH(I+3) for I = 11, 12, and 13.') 9020 FORMAT(/' Adjust I1MACH by uncommenting data statements'/ & ' appropriate for your machine.') IF (I .LT. 1 .OR. I .GT. 16) GO TO 40 I1MACH = IMACH(I) RETURN 40 WRITE(*,*) 'I1MACH(I): I =',I,' is out of bounds.' STOP ! /* C source for I1MACH -- remove the * in column 1 */ ! /* Note that some values may need changing. */ !#include !#include !#include !#include ! !long i1mach_(long *i) !{ ! switch(*i){ ! case 1: return 5; /* standard input */ ! case 2: return 6; /* standard output */ ! case 3: return 7; /* standard punch */ ! case 4: return 0; /* standard error */ ! case 5: return 32; /* bits per integer */ ! case 6: return sizeof(int); ! case 7: return 2; /* base for integers */ ! case 8: return 31; /* digits of integer base */ ! case 9: return LONG_MAX; ! case 10: return FLT_RADIX; ! case 11: return FLT_MANT_DIG; ! case 12: return FLT_MIN_EXP; ! case 13: return FLT_MAX_EXP; ! case 14: return DBL_MANT_DIG; ! case 15: return DBL_MIN_EXP; ! case 16: return DBL_MAX_EXP; ! } ! fprintf(stderr, "invalid argument: i1mach(%ld)\n", *i); ! exit(1);return 0; /* some compilers demand return values */ !} END !********************************************************************** !********************************************************************** SUBROUTINE I1MCR1(A, A1, B, C, D) !*** SPECIAL COMPUTATION FOR OLD CRAY MACHINES **** INTEGER A, A1, B, C, D A1 = 16777216*B + C A = 16777216*A1 + D END REAL FUNCTION R1MACH(I) INTEGER I ! ! SINGLE-PRECISION MACHINE CONSTANTS ! R1MACH(1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. ! R1MACH(2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. ! R1MACH(3) = B**(-T), THE SMALLEST RELATIVE SPACING. ! R1MACH(4) = B**(1-T), THE LARGEST RELATIVE SPACING. ! R1MACH(5) = LOG10(B) ! INTEGER SMALL(2) INTEGER LARGE(2) INTEGER RIGHT(2) INTEGER DIVER(2) INTEGER LOG10(2) ! needs to be (2) for AUTODOUBLE, HARRIS SLASH 6, ... INTEGER SC SAVE SMALL, LARGE, RIGHT, DIVER, LOG10, SC REAL RMACH(5) EQUIVALENCE (RMACH(1),SMALL(1)) EQUIVALENCE (RMACH(2),LARGE(1)) EQUIVALENCE (RMACH(3),RIGHT(1)) EQUIVALENCE (RMACH(4),DIVER(1)) EQUIVALENCE (RMACH(5),LOG10(1)) INTEGER J, K, L, T3E(3) DATA T3E(1) / 9777664 / DATA T3E(2) / 5323660 / DATA T3E(3) / 46980 / ! THIS VERSION ADAPTS AUTOMATICALLY TO MOST CURRENT MACHINES, ! INCLUDING AUTO-DOUBLE COMPILERS. ! TO COMPILE ON OLDER MACHINES, ADD A C IN COLUMN 1 ! ON THE NEXT LINE DATA SC/0/ ! AND REMOVE THE C FROM COLUMN 1 IN ONE OF THE SECTIONS BELOW. ! CONSTANTS FOR EVEN OLDER MACHINES CAN BE OBTAINED BY ! mail netlib@research.bell-labs.com ! send old1mach from blas ! PLEASE SEND CORRECTIONS TO dmg OR ehg@bell-labs.com. ! ! MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. ! DATA RMACH(1) / O402400000000 / ! DATA RMACH(2) / O376777777777 / ! DATA RMACH(3) / O714400000000 / ! DATA RMACH(4) / O716400000000 / ! DATA RMACH(5) / O776464202324 /, SC/987/ ! ! MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING ! 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). ! DATA SMALL(1) / 8388608 / ! DATA LARGE(1) / 2147483647 / ! DATA RIGHT(1) / 880803840 / ! DATA DIVER(1) / 889192448 / ! DATA LOG10(1) / 1067065499 /, SC/987/ ! DATA RMACH(1) / O00040000000 / ! DATA RMACH(2) / O17777777777 / ! DATA RMACH(3) / O06440000000 / ! DATA RMACH(4) / O06500000000 / ! DATA RMACH(5) / O07746420233 /, SC/987/ ! ! MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. ! DATA RMACH(1) / O000400000000 / ! DATA RMACH(2) / O377777777777 / ! DATA RMACH(3) / O146400000000 / ! DATA RMACH(4) / O147400000000 / ! DATA RMACH(5) / O177464202324 /, SC/987/ ! IF (SC .NE. 987) THEN ! *** CHECK FOR AUTODOUBLE *** SMALL(2) = 0 RMACH(1) = 1E13 IF (SMALL(2) .NE. 0) THEN ! *** AUTODOUBLED *** IF ( SMALL(1) .EQ. 1117925532 & .AND. SMALL(2) .EQ. -448790528) THEN ! *** IEEE BIG ENDIAN *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2146435071 LARGE(2) = -1 RIGHT(1) = 1017118720 RIGHT(2) = 0 DIVER(1) = 1018167296 DIVER(2) = 0 LOG10(1) = 1070810131 LOG10(2) = 1352628735 ELSE IF ( SMALL(2) .EQ. 1117925532 & .AND. SMALL(1) .EQ. -448790528) THEN ! *** IEEE LITTLE ENDIAN *** SMALL(2) = 1048576 SMALL(1) = 0 LARGE(2) = 2146435071 LARGE(1) = -1 RIGHT(2) = 1017118720 RIGHT(1) = 0 DIVER(2) = 1018167296 DIVER(1) = 0 LOG10(2) = 1070810131 LOG10(1) = 1352628735 ELSE IF ( SMALL(1) .EQ. -2065213935 & .AND. SMALL(2) .EQ. 10752) THEN ! *** VAX WITH D_FLOATING *** SMALL(1) = 128 SMALL(2) = 0 LARGE(1) = -32769 LARGE(2) = -1 RIGHT(1) = 9344 RIGHT(2) = 0 DIVER(1) = 9472 DIVER(2) = 0 LOG10(1) = 546979738 LOG10(2) = -805796613 ELSE IF ( SMALL(1) .EQ. 1267827943 & .AND. SMALL(2) .EQ. 704643072) THEN ! *** IBM MAINFRAME *** SMALL(1) = 1048576 SMALL(2) = 0 LARGE(1) = 2147483647 LARGE(2) = -1 RIGHT(1) = 856686592 RIGHT(2) = 0 DIVER(1) = 873463808 DIVER(2) = 0 LOG10(1) = 1091781651 LOG10(2) = 1352628735 ELSE WRITE(*,9010) STOP 777 END IF ELSE RMACH(1) = 1234567. IF (SMALL(1) .EQ. 1234613304) THEN ! *** IEEE *** SMALL(1) = 8388608 LARGE(1) = 2139095039 RIGHT(1) = 864026624 DIVER(1) = 872415232 LOG10(1) = 1050288283 ELSE IF (SMALL(1) .EQ. -1271379306) THEN ! *** VAX *** SMALL(1) = 128 LARGE(1) = -32769 RIGHT(1) = 13440 DIVER(1) = 13568 LOG10(1) = 547045274 ELSE IF (SMALL(1) .EQ. 1175639687) THEN ! *** IBM MAINFRAME *** SMALL(1) = 1048576 LARGE(1) = 2147483647 RIGHT(1) = 990904320 DIVER(1) = 1007681536 LOG10(1) = 1091781651 ELSE IF (SMALL(1) .EQ. 1251390520) THEN ! *** CONVEX C-1 *** SMALL(1) = 8388608 LARGE(1) = 2147483647 RIGHT(1) = 880803840 DIVER(1) = 889192448 LOG10(1) = 1067065499 ELSE DO 10 L = 1, 3 J = SMALL(1) / 10000000 K = SMALL(1) - 10000000*J IF (K .NE. T3E(L)) GO TO 20 SMALL(1) = J 10 CONTINUE ! *** CRAY T3E *** CALL I1MCRA(SMALL, K, 16, 0, 0) CALL I1MCRA(LARGE, K, 32751, 16777215, 16777215) CALL I1MCRA(RIGHT, K, 15520, 0, 0) CALL I1MCRA(DIVER, K, 15536, 0, 0) CALL I1MCRA(LOG10, K, 16339, 4461392, 10451455) GO TO 30 20 CALL I1MCRA(J, K, 16405, 9876536, 0) IF (SMALL(1) .NE. J) THEN WRITE(*,9020) STOP 777 END IF ! *** CRAY 1, XMP, 2, AND 3 *** CALL I1MCRA(SMALL(1), K, 8195, 8388608, 1) CALL I1MCRA(LARGE(1), K, 24574, 16777215, 16777214) CALL I1MCRA(RIGHT(1), K, 16338, 8388608, 0) CALL I1MCRA(DIVER(1), K, 16339, 8388608, 0) CALL I1MCRA(LOG10(1), K, 16383, 10100890, 8715216) END IF END IF 30 SC = 987 END IF ! SANITY CHECK IF (RMACH(4) .GE. 1.0) STOP 776 IF (I .LT. 1 .OR. I .GT. 5) THEN WRITE(*,*) 'R1MACH(I): I =',I,' is out of bounds.' STOP END IF R1MACH = RMACH(I) RETURN 9010 FORMAT(/' Adjust autodoubled R1MACH by getting data'/ & ' appropriate for your machine from D1MACH.') 9020 FORMAT(/' Adjust R1MACH by uncommenting data statements'/ & ' appropriate for your machine.') ! /* C source for R1MACH -- remove the * in column 1 */ !#include !#include !#include !float r1mach_(long *i) !{ ! switch(*i){ ! case 1: return FLT_MIN; ! case 2: return FLT_MAX; ! case 3: return FLT_EPSILON/FLT_RADIX; ! case 4: return FLT_EPSILON; ! case 5: return log10((double)FLT_RADIX); ! } ! fprintf(stderr, "invalid argument: r1mach(%ld)\n", *i); ! exit(1); return 0; /* else complaint of missing return value */ !} END !********************************************************************** !********************************************************************** SUBROUTINE I1MCRA(A, A1, B, C, D) !**** SPECIAL COMPUTATION FOR CRAY MACHINES **** INTEGER A, A1, B, C, D A1 = 16777216*B + C A = 16777216*A1 + D END function alog (x) ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. dimension alncs(6), center(4), alncen(5) external csevl, inits, r1mach ! ! series for aln on the interval 0. to 3.46021d-03 ! with weighted error 1.50e-16 ! log weighted error 15.82 ! significant figures required 15.65 ! decimal places required 16.21 ! data alncs( 1) / 1.3347199877973882e0 / data alncs( 2) / .000693756283284112e0 / data alncs( 3) / .000000429340390204e0 / data alncs( 4) / .000000000289338477e0 / data alncs( 5) / .000000000000205125e0 / data alncs( 6) / .000000000000000150e0 / ! data center(1) / 1.0 / data center(2) / 1.25 / data center(3) / 1.50 / data center(4) / 1.75 / ! data alncen( 1) / 0.0e0 / data alncen( 2) / .223143551314209755e0 / data alncen( 3) / .405465108108164381e0 / data alncen( 4) / .559615787935422686e0 / data alncen( 5) / .693147180559945309e0 / ! ! aln2 = alog(2.0) - 0.625 data aln2 / 0.068147180559945309e0 / data nterms / 0 / ! if (nterms.eq.0) nterms = inits (alncs, 6, 28.9*r1mach(3)) ! if (x.le.0.) call seteru ( & 29halog x is zero or negative, 29, 1, 2) ! call r9upak (x, y, n) ! 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 ! alog = 0.625*xn + (aln2*xn + alncen(ntrval) + 2.0*t + & t*t2*csevl(578.0*t2-1.0, alncs, nterms) ) ! return end !********************************************************************** !********************************************************************** function csevl (x, cs, n) ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! evaluate the n-term chebyshev series cs at x. adapted from ! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). also see fox ! and parker, chebyshev polys in numerical analysis, oxford press, p.56. ! ! input arguments -- ! x value at which the series is to be evaluated. ! cs array of n terms of a chebyshev series. in eval- ! uating cs, only half the first coef is summed. ! n number of terms in array cs. ! dimension cs(1) ! if (n.lt.1) call seteru (28hcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hcsevl number of terms gt 1000, & 31, 3, 2) if (x.lt.(-1.1) .or. x.gt.1.1) call seteru ( & 25hcsevl x outside (-1,+1), 25, 1, 1) ! 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 ! csevl = 0.5 * (b0-b2) ! return end !********************************************************************** !********************************************************************** subroutine d9gaml (xmin, xmax) ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! calculate the minimum and maximum legal bounds for x in gamma(x). ! xmin and xmax are not the only bounds, but they are the only non- ! trivial ones to calculate. ! ! output arguments -- ! xmin dble prec minimum legal value of x in gamma(x). any smaller ! value of x might result in underflow. ! xmax dble prec maximum legal value of x in gamma(x). any larger ! value of x might cause overflow. ! double precision xmin, xmax, alnbig, alnsml, xln, xold, d1mach, & dlog external d1mach, dlog ! alnsml = dlog(d1mach(1)) xmin = -alnsml do 10 i=1,10 xold = xmin xln = dlog(xmin) xmin = xmin - xmin*((xmin+0.5d0)*xln - xmin - 0.2258d0 + & alnsml) / (xmin*xln+0.5d0) if (dabs(xmin-xold).lt.0.005d0) go to 20 10 continue call seteru (27hd9gaml unable to find xmin, 27, 1, 2) ! 20 xmin = -xmin + 0.01d0 ! alnbig = dlog (d1mach(2)) xmax = alnbig do 30 i=1,10 xold = xmax xln = dlog(xmax) xmax = xmax - xmax*((xmax-0.5d0)*xln - xmax + 0.9189d0 - & alnbig) / (xmax*xln-0.5d0) if (dabs(xmax-xold).lt.0.005d0) go to 40 30 continue call seteru (27hd9gaml unable to find xmax, 27, 2, 2) ! 40 xmax = xmax - 0.01d0 xmin = dmax1 (xmin, -xmax+1.d0) ! return end !********************************************************************** !********************************************************************** double precision function d9lgmc (x) ! august 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! compute the log gamma correction factor for x .ge. 10. so that ! dlog (dgamma(x)) = dlog(dsqrt(2*pi)) + (x-.5)*dlog(x) - x + d9lgmc(x) ! double precision x, algmcs(15), xbig, xmax, dcsevl, d1mach, & dexp, dlog, dsqrt external d1mach, dcsevl, dexp, dlog, dsqrt, initds ! ! series for algm on the interval 0. to 1.00000e-02 ! with weighted error 1.28e-31 ! log weighted error 30.89 ! significant figures required 29.81 ! decimal places required 31.48 ! data algmcs( 1) / .1666389480451863247205729650822d+0 / data algmcs( 2) / -.1384948176067563840732986059135d-4 / data algmcs( 3) / .9810825646924729426157171547487d-8 / data algmcs( 4) / -.1809129475572494194263306266719d-10 / data algmcs( 5) / .6221098041892605227126015543416d-13 / data algmcs( 6) / -.3399615005417721944303330599666d-15 / data algmcs( 7) / .2683181998482698748957538846666d-17 / data algmcs( 8) / -.2868042435334643284144622399999d-19 / data algmcs( 9) / .3962837061046434803679306666666d-21 / data algmcs( 10) / -.6831888753985766870111999999999d-23 / data algmcs( 11) / .1429227355942498147573333333333d-24 / data algmcs( 12) / -.3547598158101070547199999999999d-26 / data algmcs( 13) / .1025680058010470912000000000000d-27 / data algmcs( 14) / -.3401102254316748799999999999999d-29 / data algmcs( 15) / .1276642195630062933333333333333d-30 / ! data nalgm, xbig, xmax / 0, 2*0.d0 / ! if (nalgm.ne.0) go to 10 nalgm = initds (algmcs, 15, sngl(d1mach(3)) ) xbig = 1.0d0/dsqrt(d1mach(3)) xmax = dexp (dmin1(dlog(d1mach(2)/12.d0), -dlog(12.d0*d1mach(1)))) ! 10 if (x.lt.10.d0) call seteru (23hd9lgmc x must be ge 10, 23, 1, 2) if (x.ge.xmax) go to 20 ! d9lgmc = 1.d0/(12.d0*x) if (x.lt.xbig) d9lgmc = dcsevl (2.0d0*(10.d0/x)**2-1.d0, algmcs, & nalgm) / x return ! 20 d9lgmc = 0.d0 call seteru (34hd9lgmc x so big d9lgmc underflows, 34, 2, 0) return ! end !********************************************************************** !********************************************************************** double precision function d9pak (y, n) ! december 1979 edition. w. fullerton, c3, los alamos scientific lab. ! ! pack a base 2 exponent into floating point number x. this routine is ! almost the inverse of d9upak. it is not exactly the inverse, because ! dabs(x) need not be between 0.5 and 1.0. if both d9pak and 2.d0**n ! were known to be in range we could compute ! d9pak = x * 2.0d0**n ! double precision y, aln2b, aln210, d1mach external d1mach, i1mach data nmin, nmax / 2*0 / data aln210 / 3.321928094887362347870319429489d0 / ! if (nmin.ne.0) go to 10 aln2b = 1.0d0 if (i1mach(10).ne.2) aln2b = d1mach(5)*aln210 nmin = aln2b*dble(float(i1mach(15))) nmax = aln2b*dble(float(i1mach(16))) ! 10 call d9upak (y, d9pak, ny) ! nsum = n + ny if (nsum.lt.nmin) go to 40 if (nsum.gt.nmax) call seteru ( & 31hd9pak packed number overflows, 31, 1, 2) ! if (nsum.eq.0) return if (nsum.gt.0) go to 30 ! 20 d9pak = 0.5d0*d9pak nsum = nsum + 1 if (nsum.ne.0) go to 20 return ! 30 d9pak = 2.0d0*d9pak nsum = nsum - 1 if (nsum.ne.0) go to 30 return ! 40 call seteru (32hd9pak packed number underflows, 32, 1, 0) d9pak = 0.0d0 return ! end !********************************************************************** !********************************************************************** subroutine d9upak (x, y, n) ! august 1980 portable edition. w. fullerton, los alamos scientific lab ! ! unpack floating point number x so that x = y * 2.0**n, where ! 0.5 .le. abs(y) .lt. 1.0 . ! double precision x, y, absx ! absx = dabs(x) n = 0 y = 0.0d0 if (x.eq.0.0d0) return ! 10 if (absx.ge.0.5d0) go to 20 n = n - 1 absx = absx*2.0d0 go to 10 ! 20 if (absx.lt.1.0d0) go to 30 n = n + 1 absx = absx*0.5d0 go to 20 ! 30 y = dsign (absx, x) return ! end !********************************************************************** !********************************************************************** double precision function dcsevl (x, a, n) ! ! evaluate the n-term chebyshev series a at x. adapted from ! r. broucke, algorithm 446, c.a.c.m., 16, 254 (1973). ! ! input arguments -- ! x dble prec value at which the series is to be evaluated. ! a dble prec array of n terms of a chebyshev series. in eval- ! uating a, only half the first coef is summed. ! n number of terms in array a. ! double precision a(n), x, twox, b0, b1, b2 ! if (n.lt.1) call seteru (28hdcsevl number of terms le 0, 28, 2,2) if (n.gt.1000) call seteru (31hdcsevl number of terms gt 1000, & 31, 3, 2) if (x.lt.(-1.1d0) .or. x.gt.1.1d0) call seteru ( & 25hdcsevl x outside (-1,+1), 25, 1, 1) ! twox = 2.0d0*x b1 = 0.d0 b0 = 0.d0 do 10 i=1,n b2 = b1 b1 = b0 ni = n - i + 1 b0 = twox*b1 - b2 + a(ni) 10 continue ! dcsevl = 0.5d0 * (b0-b2) ! return end !********************************************************************** !********************************************************************** double precision function dexp (x) ! may 1980 edition. w. fullerton, c3, los alamos scientific lab. double precision x, expcs(14), twon16(17), aln216, f, xint, xmax,& xmin, y, d1mach, dint, d9pak, dcsevl, dlog external d1mach, d9pak, dcsevl, dint, dlog, initds ! ! series for exp on the interval -1.00000e+00 to 1.00000e+00 ! with weighted error 2.30e-34 ! log weighted error 33.64 ! significant figures required 32.28 ! decimal places required 34.21 ! data expcs( 1) / .866569493314985712733404647266231d-1 / data expcs( 2) / .938494869299839561896336579701203d-3 / data expcs( 3) / .677603970998168264074353014653601d-5 / data expcs( 4) / .366931200393805927801891250687610d-7 / data expcs( 5) / .158959053617461844641928517821508d-9 / data expcs( 6) / .573859878630206601252990815262106d-12 / data expcs( 7) / .177574448591421511802306980226000d-14 / data expcs( 8) / .480799166842372422675950244533333d-17 / data expcs( 9) / .115716376881828572809260000000000d-19 / data expcs( 10) / .250650610255497719932458666666666d-22 / data expcs( 11) / .493571708140495828480000000000000d-25 / data expcs( 12) / .890929572740634240000000000000000d-28 / data expcs( 13) / .148448062907997866666666666666666d-30 / data expcs( 14) / .229678916630186666666666666666666d-33 / ! ! twon16(i) is 2.0**((i-1)/16) - 1.0 data twon16( 1) / 0.0d0 / data twon16( 2) / .44273782427413840321966478739929d-1 / data twon16( 3) / .90507732665257659207010655760707d-1 / data twon16( 4) / .13878863475669165370383028384151d0 / data twon16( 5) / .18920711500272106671749997056047d0 / data twon16( 6) / .24185781207348404859367746872659d0 / data twon16( 7) / .29683955465100966593375411779245d0 / data twon16( 8) / .35425554693689272829801474014070d0 / data twon16( 9) / .41421356237309504880168872420969d0 / data twon16( 10) / .47682614593949931138690748037404d0 / data twon16( 11) / .54221082540794082361229186209073d0 / data twon16( 12) / .61049033194925430817952066735740d0 / data twon16( 13) / .68179283050742908606225095246642d0 / data twon16( 14) / .75625216037329948311216061937531d0 / data twon16( 15) / .83400808640934246348708318958828d0 / data twon16( 16) / .91520656139714729387261127029583d0 / data twon16( 17) / 1.d0 / ! ! aln216 is 16.0/alog(2.) - 23.0 data aln216 / .83120654223414517758794896030274d-1 / data nterms, xmin, xmax /0, 2*0.d0 / ! if (nterms.ne.0) go to 10 nterms = initds (expcs, 14, 0.1*sngl(d1mach(3))) xmin = dlog (d1mach(1)) + .01d0 xmax = dlog (d1mach(2)) - 0.001d0 ! 10 if (x.lt.xmin) go to 20 if (x.gt.xmax) call seteru ( & 31hdexp x so big dexp overflows, 31, 2, 2) ! xint = dint (x) y = x - xint ! y = 23.d0*y + x*aln216 n = y f = y - dble(float(n)) n = 23.d0*xint + dble(float(n)) n16 = n/16 if (n.lt.0) n16 = n16 - 1 ndx = n - 16*n16 + 1 ! dexp = 1.0d0 + (twon16(ndx) + f*(1.0d0 + twon16(ndx)) * & dcsevl (f, expcs, nterms) ) ! dexp = d9pak (dexp, n16) return ! 20 call seteru (34hdexp x so small dexp underflows, 34, 1, 0) dexp = 0.d0 return ! end !********************************************************************** !********************************************************************** double precision function dgamma (x) ! jan 1984 edition. w. fullerton, c3, los alamos scientific lab. ! jan 1994 wpp@ips.id.ethz.ch, ehg@research.att.com declare xsml double precision x, gamcs(42), dxrel, pi, sinpiy, sq2pil, xmax, & xmin, y, d9lgmc, dcsevl, d1mach, dexp, dint, dlog, & dsin, dsqrt, xsml external d1mach, d9lgmc, dcsevl, dexp, dint, dlog, dsin, dsqrt, & initds ! ! series for gam on the interval 0. to 1.00000e+00 ! with weighted error 5.79e-32 ! log weighted error 31.24 ! significant figures required 30.00 ! decimal places required 32.05 ! data gamcs( 1) / .8571195590989331421920062399942d-2 / data gamcs( 2) / .4415381324841006757191315771652d-2 / data gamcs( 3) / .5685043681599363378632664588789d-1 / data gamcs( 4) / -.4219835396418560501012500186624d-2 / data gamcs( 5) / .1326808181212460220584006796352d-2 / data gamcs( 6) / -.1893024529798880432523947023886d-3 / data gamcs( 7) / .3606925327441245256578082217225d-4 / data gamcs( 8) / -.6056761904460864218485548290365d-5 / data gamcs( 9) / .1055829546302283344731823509093d-5 / data gamcs( 10) / -.1811967365542384048291855891166d-6 / data gamcs( 11) / .3117724964715322277790254593169d-7 / data gamcs( 12) / -.5354219639019687140874081024347d-8 / data gamcs( 13) / .9193275519859588946887786825940d-9 / data gamcs( 14) / -.1577941280288339761767423273953d-9 / data gamcs( 15) / .2707980622934954543266540433089d-10 / data gamcs( 16) / -.4646818653825730144081661058933d-11 / data gamcs( 17) / .7973350192007419656460767175359d-12 / data gamcs( 18) / -.1368078209830916025799499172309d-12 / data gamcs( 19) / .2347319486563800657233471771688d-13 / data gamcs( 20) / -.4027432614949066932766570534699d-14 / data gamcs( 21) / .6910051747372100912138336975257d-15 / data gamcs( 22) / -.1185584500221992907052387126192d-15 / data gamcs( 23) / .2034148542496373955201026051932d-16 / data gamcs( 24) / -.3490054341717405849274012949108d-17 / data gamcs( 25) / .5987993856485305567135051066026d-18 / data gamcs( 26) / -.1027378057872228074490069778431d-18 / data gamcs( 27) / .1762702816060529824942759660748d-19 / data gamcs( 28) / -.3024320653735306260958772112042d-20 / data gamcs( 29) / .5188914660218397839717833550506d-21 / data gamcs( 30) / -.8902770842456576692449251601066d-22 / data gamcs( 31) / .1527474068493342602274596891306d-22 / data gamcs( 32) / -.2620731256187362900257328332799d-23 / data gamcs( 33) / .4496464047830538670331046570666d-24 / data gamcs( 34) / -.7714712731336877911703901525333d-25 / data gamcs( 35) / .1323635453126044036486572714666d-25 / data gamcs( 36) / -.2270999412942928816702313813333d-26 / data gamcs( 37) / .3896418998003991449320816639999d-27 / data gamcs( 38) / -.6685198115125953327792127999999d-28 / data gamcs( 39) / .1146998663140024384347613866666d-28 / data gamcs( 40) / -.1967938586345134677295103999999d-29 / data gamcs( 41) / .3376448816585338090334890666666d-30 / data gamcs( 42) / -.5793070335782135784625493333333d-31 / ! data pi / 3.14159265358979323846264338327950d0 / ! sq2pil is 0.5*alog(2*pi) = alog(sqrt(2*pi)) data sq2pil / 0.91893853320467274178032973640562d0 / data ngam, xmin, xmax, xsml, dxrel / 0, 4*0.d0 / ! if (ngam.ne.0) go to 10 ngam = initds (gamcs, 42, 0.1*sngl(d1mach(3)) ) ! call d9gaml (xmin, xmax) xsml = dexp (dmax1 (dlog(d1mach(1)), -dlog(d1mach(2)))+0.01d0) dxrel = dsqrt (d1mach(4)) ! 10 y = dabs(x) if (y.gt.10.d0) go to 50 ! ! compute gamma(x) for -xbnd .le. x .le. xbnd. reduce interval and find ! gamma(1+y) for 0.0 .le. y .lt. 1.0 first of all. ! n = x if (x.lt.0.d0) n = n - 1 y = x - dble(float(n)) n = n - 1 dgamma = 0.9375d0 + dcsevl (2.d0*y-1.d0, gamcs, ngam) if (n.eq.0) return ! if (n.gt.0) go to 30 ! ! compute gamma(x) for x .lt. 1.0 ! n = -n if (x.eq.0.d0) call seteru (14hdgamma x is 0, 14, 4, 2) if (x.lt.0.0d0 .and. x+dble(float(n-2)).eq.0.d0) call seteru ( & 31hdgamma x is a negative integer, 31, 4, 2) if (x.lt.(-0.5d0) .and. dabs((x-dint(x-0.5d0))/x).lt.dxrel) & call seteru ('68hdgamma answer lt half precision because x' // & 'too near negative integer', 68, 1, 1) if (y.lt.xsml) call seteru ( & 54hdgamma x is so close to 0.0 that the result overflows, & 54, 5, 2) ! do 20 i=1,n dgamma = dgamma/(x+dble(float(i-1)) ) 20 continue return ! ! gamma(x) for x .ge. 2.0 and x .le. 10.0 ! 30 do 40 i=1,n dgamma = (y+dble(float(i))) * dgamma 40 continue return ! ! gamma(x) for dabs(x) .gt. 10.0. recall y = dabs(x). ! 50 if (x.gt.xmax) call seteru (32hdgamma x so big gamma overflows, & 32, 3, 2) ! dgamma = 0.d0 if (x.lt.xmin) call seteru (35hdgamma x so small gamma underflows, & 35, 2, 0) if (x.lt.xmin) return ! dgamma = dexp ((y-0.5d0)*dlog(y) - y + sq2pil + d9lgmc(y) ) if (x.gt.0.d0) return ! if (dabs((x-dint(x-0.5d0))/x).lt.dxrel) call seteru ( & '61hdgamma answer lt half precision, x too near negative' // & 'integer' , 61, 1, 1) ! sinpiy = dsin (pi*y) if (sinpiy.eq.0.d0) call seteru ( & 31hdgamma x is a negative integer, 31, 4, 2) ! dgamma = -pi/(y*sinpiy*dgamma) ! return end !********************************************************************** !********************************************************************** double precision function dint (x) ! december 1983 edition. w. fullerton, c3, los alamos scientific lab. ! ! dint is the double precision equivalent of aint. this portable ! version is quite efficient when the argument is reasonably small (a ! common case), and so no faster machine-dependent version is needed. ! double precision x, xscl, scale, xbig, xmax, part, d1mach, & dlog external d1mach, dlog, i1mach, r1mach data npart, scale, xbig, xmax / 0, 3*0.0d0 / ! if (npart.ne.0) go to 10 ibase = i1mach(10) xmax = 1.0d0/d1mach (4) xbig = amin1 (float (i1mach(9)), 1.0/r1mach(4)) scale = ibase**int(dlog(xbig)/dlog(dble(float(ibase)))-0.5d0) npart = dlog(xmax)/dlog(scale) + 1.0d0 ! 10 if (x.lt.(-xbig) .or. x.gt.xbig) go to 20 ! dint = int(sngl(x)) return ! 20 xscl = dabs(x) if (xscl.gt.xmax) go to 50 ! do 30 i=1,npart xscl = xscl/scale 30 continue ! dint = 0.0d0 do 40 i=1,npart xscl = xscl*scale ipart = xscl part = ipart xscl = xscl - part dint = dint*scale + part 40 continue ! if (x.lt.0.0d0) dint = -dint return ! 50 call seteru ('68hdint dabs(x) may be too big to be represented' // & 'as an exact integer', 68, 1, 1) dint = x return ! end !********************************************************************** !********************************************************************** double precision function dlog (x) ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, alncs(11), center(4), alncen(5), aln2, y, t, & t2, xn, dcsevl, d1mach external d1mach, dcsevl, initds ! ! series for aln on the interval 0. to 3.46021e-03 ! with weighted error 4.15e-32 ! log weighted error 31.38 ! significant figures required 31.21 ! decimal places required 31.90 ! data alncs( 1) / .13347199877973881561689386047187d+1 / data alncs( 2) / .69375628328411286281372438354225d-3 / data alncs( 3) / .42934039020450834506559210803662d-6 / data alncs( 4) / .28933847795432594580466440387587d-9 / data alncs( 5) / .20512517530340580901741813447726d-12 / data alncs( 6) / .15039717055497386574615153319999d-15 / data alncs( 7) / .11294540695636464284521613333333d-18 / data alncs( 8) / .86355788671171868881946666666666d-22 / data alncs( 9) / .66952990534350370613333333333333d-25 / data alncs( 10) / .52491557448151466666666666666666d-28 / data alncs( 11) / .41530540680362666666666666666666d-31 / ! data center(1) / 1.0d0 / data center(2) / 1.25d0 / data center(3) / 1.50d0 / data center(4) / 1.75d0 / ! data alncen( 1) / 0.0d0 / data alncen( 2) / .22314355131420975576629509030983d0 / data alncen( 3) / .40546510810816438197801311546434d0 / data alncen( 4) / .55961578793542268627088850052682d0 / data alncen( 5) / .69314718055994530941723212145817d0 / ! ! aln2 = alog(2.0) - 0.625 data aln2 / 0.06814718055994530941723212145818d0 / data nterms / 0 / ! if (nterms.eq.0) nterms = initds (alncs, 11, 28.9*sngl(d1mach(3))) ! if (x.le.0.d0) call seteru ( & 29hdlog x is zero or negative, 29, 1, 2) ! call d9upak (x, y, n) ! xn = n - 1 y = 2.0d0*y ntrval = 4.0d0*y - 2.5d0 ! if (ntrval.eq.5) t = ((y-1.0d0)-1.0d0) / (y+2.0d0) if (ntrval.lt.5) t = (y-center(ntrval)) / (y+center(ntrval)) t2 = t*t dlog = 0.625d0*xn + (aln2*xn + alncen(ntrval) + 2.0d0*t + & t*t2*dcsevl(578.d0*t2-1.0d0, alncs, nterms) ) ! return end !********************************************************************** !********************************************************************** double precision function dsin (x) ! august 1980 edition. w. fullerton, los alamos scientific lab. ! ! this routine is based on the algorithm of cody and waite in ! argonne tm-321, software manual working note number 1 ! double precision x, sincs(15), pihi, pilo, pirec, pi2rec, xsml, & xwarn, xmax, y, xn, sgn, f, dint, dcsevl, d1mach, dsqrt external d1mach, dcsevl, dint, dsqrt, initds ! ! series for sin on the interval 0.00000e+00 to 2.46740e+00 ! with weighted error 2.56e-34 ! log weighted error 33.59 ! significant figures required 33.01 ! decimal places required 34.18 ! data sincs( 1) / -0.374991154955873175839919279977323464d0/ data sincs( 2) / -0.181603155237250201863830316158004754d0/ data sincs( 3) / 0.005804709274598633559427341722857921d0/ data sincs( 4) / -0.000086954311779340757113212316353178d0/ data sincs( 5) / 0.000000754370148088851481006839927030d0/ data sincs( 6) / -0.000000004267129665055961107126829906d0/ data sincs( 7) / 0.000000000016980422945488168181824792d0/ data sincs( 8) / -0.000000000000050120578889961870929524d0/ data sincs( 9) / 0.000000000000000114101026680010675628d0/ data sincs( 10) / -0.000000000000000000206437504424783134d0/ data sincs( 11) / 0.000000000000000000000303969595918706d0/ data sincs( 12) / -0.000000000000000000000000371357734157d0/ data sincs( 13) / 0.000000000000000000000000000382486123d0/ data sincs( 14) / -0.000000000000000000000000000000336623d0/ data sincs( 15) / 0.000000000000000000000000000000000256d0/ ! ! pihi + pilo = pi. pihi is exactly representable on all machines ! with at least 8 bits of precision. whether it is exactly ! represented depends on the compiler. this routine is more ! accurate if it is exactly represented. data pihi / 3.140625d0 / data pilo / 9.6765358979323846264338327950288d-4/ data pirec / 0.31830988618379067153776752674503d0 / data pi2rec / 0.63661977236758134307553505349006d0 / data ntsn, xsml, xwarn, xmax / 0, 3*0.0d0 / ! if (ntsn.ne.0) go to 10 ntsn = initds (sincs, 15, 0.1*sngl(d1mach(3))) ! xsml = dsqrt (2.0d0*d1mach(3)) xmax = 1.0d0/d1mach(4) xwarn = dsqrt (xmax) ! 10 y = dabs (x) if (y.gt.xmax) call seteru ( & 42hdsin no precision because abs(x) is big, 42, 2, 2) if (y.gt.xwarn) call seteru ( & 54hdsin answer lt half precision because abs(x) is big, & 54, 1, 1) ! dsin = x if (y.lt.xsml) return ! xn = dint (y*pirec+0.5d0) n2 = dmod (xn, 2.0d0) + 0.5d0 sgn = x if (n2.ne.0) sgn = -sgn f = (y-xn*pihi) - xn*pilo ! dsin = f + f*dcsevl(2.0d0*(f*pi2rec)**2-1.0d0, sincs, ntsn) if (sgn.lt.0.0d0) dsin = -dsin if (dabs(dsin).gt.1.0d0) dsin = dsign (1.0d0, dsin) ! return end !********************************************************************** !********************************************************************** double precision function dsqrt (x) ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. double precision x, sqrt2(3), y, d9pak, d1mach external alog, d1mach, d9pak data sqrt2(1) / 0.70710678118654752440084436210485d0 / data sqrt2(2) / 1.0d0 / data sqrt2(3) / 1.41421356237309504880168872420970d0 / ! data niter / 0 / ! if (niter.eq.0) niter = 1.443*alog(-0.104*alog(0.1* & sngl(d1mach(3)) )) + 1.0 ! if (x.le.0.d0) go to 20 ! call d9upak (x, y, n) ixpnt = n/2 irem = n - 2*ixpnt + 2 ! ! the approximation below has accuracy of 4.16 digits. z = y dsqrt = .261599e0 + z*(1.114292e0 + z*(-.516888e0 + z*.141067e0)) ! do 10 iter=1,niter dsqrt = dsqrt + 0.5d0*(y - dsqrt*dsqrt) / dsqrt 10 continue ! dsqrt = d9pak (sqrt2(irem)*dsqrt, ixpnt) return ! 20 if (x.lt.0.d0) call seteru (21hdsqrt x is negative, 21, 1, 1) dsqrt = 0.0d0 return ! end !********************************************************************** !********************************************************************** subroutine e9rint(messg,nw,nerr,save) ! ! this routine stores the current error message or prints the old one, ! if any, depending on whether or not save = .true. . ! integer messg(nw) logical save external i1mach, i8save ! ! messgp stores at least the first 72 characters of the previous ! message. its length is machine dependent and must be at least ! ! 1 + 71/(the number of characters stored per integer word). ! integer messgp(36),fmt(14),ccplus ! ! start with no previous message. ! data messgp(1)/1h1/, nwp/0/, nerrp/0/ ! ! set up the format for printing the error message. ! the format is simply (a1,14x,72axx) where xx=i1mach(6) is the ! number of characters stored per integer word. ! data ccplus / 1h+ / ! data fmt( 1) / 1h( / data fmt( 2) / 1ha / data fmt( 3) / 1h1 / data fmt( 4) / 1h, / data fmt( 5) / 1h1 / data fmt( 6) / 1h4 / data fmt( 7) / 1hx / data fmt( 8) / 1h, / data fmt( 9) / 1h7 / data fmt(10) / 1h2 / data fmt(11) / 1ha / data fmt(12) / 1hx / data fmt(13) / 1hx / data fmt(14) / 1h) / ! if (.not.save) go to 20 ! ! save the message. ! nwp=nw nerrp=nerr do 10 i=1,nw 10 messgp(i)=messg(i) ! go to 30 ! 20 if (i8save(1,0,.false.).eq.0) go to 30 ! ! print the message. ! iwunit=i1mach(4) write(iwunit,9000) nerrp 9000 format(7h error ,i4,4h in ) ! call s88fmt(2,i1mach(6),fmt(12)) write(iwunit,fmt) ccplus,(messgp(i),i=1,nwp) ! 30 return ! end !********************************************************************** !********************************************************************** subroutine eprint ! ! this subroutine prints the last error message, if any. ! integer messg(1) ! call e9rint(messg,1,1,.false.) return ! end !********************************************************************** !********************************************************************** integer function i8save(isw,ivalue,set) ! ! if (isw = 1) i8save returns the current error number and ! sets it to ivalue if set = .true. . ! ! if (isw = 2) i8save returns the current recovery switch and ! sets it to ivalue if set = .true. . ! logical set ! integer iparam(2) ! iparam(1) is the error number and iparam(2) is the recovery switch. ! ! start execution error free and with recovery turned off. ! data iparam(1) /0/, iparam(2) /2/ ! i8save=iparam(isw) if (set) iparam(isw)=ivalue ! return ! end !********************************************************************** !********************************************************************** function initds (dos, nos, eta) ! june 1977 edition. w. fullerton, c3, los alamos scientific lab. ! ! initialize the double precision orthogonal series dos so that initds ! is the number of terms needed to insure the error is no larger than ! eta. ordinarily eta will be chosen to be one-tenth machine precision. ! ! input arguments -- ! dos dble prec array of nos coefficients in an orthogonal series. ! nos number of coefficients in dos. ! eta requested accuracy of series. ! double precision dos(nos) ! if (nos.lt.1) call seteru ( & 35hinitds number of coefficients lt 1, 35, 2, 2) ! err = 0. do 10 ii=1,nos i = nos + 1 - ii err = err + abs(sngl(dos(i))) if (err.gt.eta) go to 20 10 continue ! 20 if (i.eq.nos) call seteru (28hinitds eta may be too small, 28, & 1, 2) initds = i ! return end !********************************************************************** !********************************************************************** function inits (os, nos, eta) ! april 1977 version. w. fullerton, c3, los alamos scientific lab. ! ! initialize the orthogonal series so that inits is the number of terms ! needed to insure the error is no larger than eta. ordinarily, eta ! will be chosen to be one-tenth machine precision. ! ! input arguments -- ! os array of nos coefficients in an orthogonal series. ! nos number of coefficients in os. ! eta requested accuracy of series. ! dimension os(nos) ! if (nos.lt.1) call seteru ( & 35hinits number of coefficients lt 1, 35, 2, 2) ! 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 ! 20 if (i.eq.nos) call seteru (28hinits eta may be too small, 28, & 1, 2) inits = i ! return end !********************************************************************** !********************************************************************** subroutine r9upak (x, y, n) ! august 1980 portable edition. w. fullerton, los alamos scientific lab ! ! unpack floating point number x so that x = y * 2.0**n, where ! 0.5 .le. abs(y) .lt. 1.0 . ! absx = abs(x) n = 0 y = 0.0 if (x.eq.0.0) return ! 10 if (absx.ge.0.5) go to 20 n = n - 1 absx = absx*2.0 go to 10 ! 20 if (absx.lt.1.0) go to 30 n = n + 1 absx = absx*0.5 go to 20 ! 30 y = sign (absx, x) return ! end !********************************************************************** !********************************************************************** subroutine s88fmt( n, w, ifmt ) ! ! s88fmt replaces ifmt(1), ... , ifmt(n) with ! the characters corresponding to the n least significant ! digits of w. ! integer n,w,ifmt(n) ! integer nt,wt,digits(10) ! data digits( 1) / 1h0 / data digits( 2) / 1h1 / data digits( 3) / 1h2 / data digits( 4) / 1h3 / data digits( 5) / 1h4 / data digits( 6) / 1h5 / data digits( 7) / 1h6 / data digits( 8) / 1h7 / data digits( 9) / 1h8 / data digits(10) / 1h9 / ! nt = n wt = w ! 10 if (nt .le. 0) return idigit = mod( wt, 10 ) ifmt(nt) = digits(idigit+1) wt = wt/10 nt = nt - 1 go to 10 ! end !********************************************************************** !********************************************************************** subroutine seterr (messg, nmessg, nerr, iopt) ! ! this version modified by w. fullerton to dump if iopt = 1 and ! not recovering. ! seterr sets lerror = nerr, optionally prints the message and dumps ! according to the following rules... ! ! if iopt = 1 and recovering - just remember the error. ! if iopt = 1 and not recovering - print, dump and stop. ! if iopt = 2 - print, dump and stop. ! ! input ! ! messg - the error message. ! nmessg - the length of the message, in characters. ! nerr - the error number. must have nerr non-zero. ! iopt - the option. must have iopt=1 or 2. ! ! error states - ! ! 1 - message length not positive. ! 2 - cannot have nerr=0. ! 3 - an unrecovered error followed by another error. ! 4 - bad value for iopt. ! ! only the first 72 characters of the message are printed. ! ! the error handler calls a subroutine named fdump to produce a ! symbolic dump. to complete the package, a dummy version of fdump ! is supplied, but it should be replaced by a locally written version ! which at least gives a trace-back. ! integer messg(1) external i1mach, i8save ! ! the unit for error messages. ! iwunit=i1mach(4) ! if (nmessg.ge.1) go to 10 ! ! a message of non-positive length is fatal. ! write(iwunit,9000) 9000 format(52h1error 1 in seterr - message length not positive.) go to 60 ! ! nw is the number of words the message occupies. ! 10 nw=(min0(nmessg,72)-1)/i1mach(6)+1 ! if (nerr.ne.0) go to 20 ! ! cannot turn the error state off using seterr. ! write(iwunit,9001) 9001 format(42h1error 2 in seterr - cannot have nerr=0// & 34h the current error message follows///) call e9rint(messg,nw,nerr,.true.) itemp=i8save(1,1,.true.) go to 50 ! ! set lerror and test for a previous unrecovered error. ! 20 if (i8save(1,nerr,.true.).eq.0) go to 30 ! write(iwunit,9002) 9002 format(23h1error 3 in seterr -, & 48h an unrecovered error followed by another error.// & 48h the previous and current error messages follow.///) call eprint call e9rint(messg,nw,nerr,.true.) go to 50 ! ! save this message in case it is not recovered from properly. ! 30 call e9rint(messg,nw,nerr,.true.) ! if (iopt.eq.1 .or. iopt.eq.2) go to 40 ! ! must have iopt = 1 or 2. ! write(iwunit,9003) 9003 format(42h1error 4 in seterr - bad value for iopt// & 34h the current error message follows///) go to 50 ! ! test for recovery. ! 40 if (iopt.eq.2) go to 50 ! if (i8save(2,0,.false.).eq.1) return ! ! call eprint ! stop ! 50 call eprint 60 call fdump stop ! end !********************************************************************** !********************************************************************** subroutine seteru (messg, nmessg, nerr, iopt) common /cseter/ iunflo integer messg(1) data iunflo / 0 / ! if (iopt.ne.0) call seterr (messg, nmessg, nerr, iopt) if (iopt.ne.0) return ! if (iunflo.le.0) return call seterr (messg, nmessg, nerr, 1) ! return end !********************************************************************** !********************************************************************** !DECK FDUMP SUBROUTINE FDUMP !***BEGIN PROLOGUE FDUMP !***PURPOSE Symbolic dump (should be locally written). !***LIBRARY SLATEC (XERROR) !***CATEGORY R3 !***TYPE ALL (FDUMP-A) !***KEYWORDS ERROR, XERMSG !***AUTHOR Jones, R. E., (SNLA) !***DESCRIPTION ! ! ***Note*** Machine Dependent Routine ! FDUMP is intended to be replaced by a locally written ! version which produces a symbolic dump. Failing this, ! it should be replaced by a version which prints the ! subprogram nesting list. Note that this dump must be ! printed on each of up to five files, as indicated by the ! XGETUA routine. See XSETUA and XGETUA for details. ! ! Written by Ron Jones, with SLATEC Common Math Library Subcommittee ! !***REFERENCES (NONE) !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 790801 DATE WRITTEN ! 861211 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) !***END PROLOGUE FDUMP !***FIRST EXECUTABLE STATEMENT FDUMP RETURN END !*********************************************************************** !*********************************************************************** SUBROUTINE INVERT(N,A,AINV,INV_IO) !INPUT : N,A,INV_IO !OUTPUT: AINV !THIS PROGRAM COMPUTES THE INVERSE OF THE NxN MATRIX A IF IT IS NOT ! SINGULAR !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY INVERT*********************************** ! MATMUL !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY INVERT************************************ ! DGEFA,DGECO,DGEDI !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER N,INV_IO DOUBLE PRECISION, DIMENSION(N,N) :: & A !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N,N) :: & AINV !INTERNAL VARIABLES INTEGER I,J,INFO,JOB INTEGER, DIMENSION(N) :: & IPVT DOUBLE PRECISION RCOND DOUBLE PRECISION DET(2) DOUBLE PRECISION, DIMENSION(N) :: & Z,WORK DOUBLE PRECISION, DIMENSION(N,N) :: & ATEMP,PROD !COMPUTE ONLY THE L*U FACTORIZATION OF MATRIX A FOR SUBROUTINE DGEDI ! IF NOT DOING DIAGNOSTIC TESTING (I.E. NORMAL OPERATION) IF (INV_IO == 0) THEN ATEMP = A CALL DGEFA(ATEMP,N,N,IPVT,INFO) !COMPUTE THE L*U FACTORIZATION OF MATRIX A FOR SUBROUTINE !DGEDI AND ITS RECIPROCAL CONDITION NUMBER FOR DIAGNOSTIC !TESTING PURPOSES ELSE ATEMP = A CALL DGECO(ATEMP,N,N,IPVT,RCOND,Z) PRINT* PRINT*,'THE RECIPROCAL CONDITION NUMBER OF THE MATRIX IS:' PRINT*, RCOND END IF !COMPUTE THE INVERSE OF MATRIX A USING OUTPUT OF SUBROUTINE DGEFA OR ! DGECO JOB = 1 CALL DGEDI(ATEMP,N,N,IPVT,DET,WORK,JOB) AINV = ATEMP IF (INV_IO == 1) THEN PRINT* PRINT*,'THE MATRIX INVERSE LOOKS LIKE:' WRITE(*,10) ((AINV(I,J),J=1,N),I=1,N) 10 FORMAT(8(1X,F14.8)) PROD = MATMUL(A,AINV) PRINT* PRINT*,'THE PRODUCT OF THE ORIGINAL MATRIX' PRINT*,'AND ITS INVERSE LOOKS LIKE (SHOULD BE IDENTITY):' WRITE(*,10) ((PROD(I,J),J=1,N),I=1,N) END IF RETURN END !*********************************************************************** !*********************************************************************** !DECK DASUM DOUBLE PRECISION FUNCTION DASUM (N, DX, INCX)! !***BEGIN PROLOGUE DASUM !***PURPOSE Compute the sum of the magnitudes of the elements of a ! vector. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A3A !***TYPE DOUBLE PRECISION (SASUM-S, DASUM-D, SCASUM-C) !***KEYWORDS BLAS, LINEAR ALGEBRA, SUM OF MAGNITUDES OF A VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! ! --Output-- ! DASUM double precision result (zero if N .LE. 0) ! ! Returns sum of magnitudes of double precision DX. ! DASUM = sum from 0 to N-1 of ABS(DX(IX+I*INCX)), ! where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900821 Modified to correct problem with a negative increment. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DASUM IMPLICIT NONE DOUBLE PRECISION DX(*) INTEGER I, INCX, IX, M, MP1, N !***FIRST EXECUTABLE STATEMENT DASUM DASUM = 0.0D0 IF (N .LE. 0) RETURN ! IF (INCX .EQ. 1) GOTO 20 ! ! Code for increment not equal to 1. ! IX = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 DO 10 I = 1,N DASUM = DASUM + ABS(DX(IX)) IX = IX + INCX 10 CONTINUE RETURN ! ! Code for increment equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 6. ! 20 M = MOD(N,6) IF (M .EQ. 0) GOTO 40 DO 30 I = 1,M DASUM = DASUM + ABS(DX(I)) 30 CONTINUE IF (N .LT. 6) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,6 DASUM = DASUM + ABS(DX(I)) + ABS(DX(I+1)) + ABS(DX(I+2)) + & ABS(DX(I+3)) + ABS(DX(I+4)) + ABS(DX(I+5)) 50 CONTINUE RETURN END !************************************************************************** !*************************************************************************** !DECK DAXPY SUBROUTINE DAXPY (N, DA, DX, INCX, DY, INCY) !***BEGIN PROLOGUE DAXPY !***PURPOSE Compute a constant times a vector plus a vector. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A7 !***TYPE DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C) !***KEYWORDS BLAS, LINEAR ALGEBRA, TRIAD, VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DA double precision scalar multiplier ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! DY double precision vector with N elements ! INCY storage spacing between elements of DY ! ! --Output-- ! DY double precision result (unchanged if N .LE. 0) ! ! Overwrite double precision DY with double precision DA*DX + DY. ! For I = 0 to N-1, replace DY(LY+I*INCY) with DA*DX(LX+I*INCX) + ! DY(LY+I*INCY), ! where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is ! defined in a similar way using INCY. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920310 Corrected definition of LX in DESCRIPTION. (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DAXPY IMPLICIT NONE INTEGER N,INCX,INCY,IX,IY,I,M,MP1,NS DOUBLE PRECISION DX(*), DY(*), DA !***FIRST EXECUTABLE STATEMENT DAXPY IF (N.LE.0 .OR. DA.EQ.0.0D0) RETURN IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60 ! ! Code for unequal or nonpositive increments. ! 5 IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN ! ! Code for both increments equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 4. ! 20 M = MOD(N,4) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF (N .LT. 4) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I+1) = DY(I+1) + DA*DX(I+1) DY(I+2) = DY(I+2) + DA*DX(I+2) DY(I+3) = DY(I+3) + DA*DX(I+3) 50 CONTINUE RETURN ! ! Code for equal, positive, non-unit increments. ! 60 NS = N*INCX DO 70 I = 1,NS,INCX DY(I) = DA*DX(I) + DY(I) 70 CONTINUE RETURN END !*************************************************************************** !*************************************************************************** !DECK DDOT DOUBLE PRECISION FUNCTION DDOT (N, DX, INCX, DY, INCY) !***BEGIN PROLOGUE DDOT !***PURPOSE Compute the inner product of two vectors. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A4 !***TYPE DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C) !***KEYWORDS BLAS, INNER PRODUCT, LINEAR ALGEBRA, VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! DY double precision vector with N elements ! INCY storage spacing between elements of DY ! ! --Output-- ! DDOT double precision dot product (zero if N .LE. 0) ! ! Returns the dot product of double precision DX and DY. ! DDOT = sum for I = 0 to N-1 of DX(LX+I*INCX) * DY(LY+I*INCY), ! where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is ! defined in a similar way using INCY. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920310 Corrected definition of LX in DESCRIPTION. (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DDOT IMPLICIT NONE INTEGER N,INCX,INCY,IX,IY,I,M,MP1,NS DOUBLE PRECISION DX(*), DY(*) !***FIRST EXECUTABLE STATEMENT DDOT DDOT = 0.0D0 IF (N .LE. 0) RETURN IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60 ! ! Code for unequal or nonpositive increments. ! 5 IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DDOT = DDOT + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN ! ! Code for both increments equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 5. ! 20 M = MOD(N,5) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DDOT = DDOT + DX(I)*DY(I) 30 CONTINUE IF (N .LT. 5) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + DX(I+2)*DY(I+2) + & DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4) 50 CONTINUE RETURN ! ! Code for equal, positive, non-unit increments. ! 60 NS = N*INCX DO 70 I = 1,NS,INCX DDOT = DDOT + DX(I)*DY(I) 70 CONTINUE RETURN END !*************************************************************************** !*************************************************************************** !DECK DGECO SUBROUTINE DGECO (A, LDA, N, IPVT, RCOND, Z) !***BEGIN PROLOGUE DGECO !***PURPOSE Factor a matrix using Gaussian elimination and estimate ! the condition number of the matrix. !***LIBRARY SLATEC (LINPACK) !***CATEGORY D2A1 !***TYPE DOUBLE PRECISION (SGECO-S, DGECO-D, CGECO-C) !***KEYWORDS CONDITION NUMBER, GENERAL MATRIX, LINEAR ALGEBRA, LINPACK, ! MATRIX FACTORIZATION !***AUTHOR Moler, C. B., (U. of New Mexico) !***DESCRIPTION ! ! DGECO factors a double precision matrix by Gaussian elimination ! and estimates the condition of the matrix. ! ! If RCOND is not needed, DGEFA is slightly faster. ! To solve A*X = B , follow DGECO by DGESL. ! To compute INVERSE(A)*C , follow DGECO by DGESL. ! To compute DETERMINANT(A) , follow DGECO by DGEDI. ! To compute INVERSE(A) , follow DGECO by DGEDI. ! ! On Entry ! ! A DOUBLE PRECISION(LDA, N) ! the matrix to be factored. ! ! LDA INTEGER ! the leading dimension of the array A . ! ! N INTEGER ! the order of the matrix A . ! ! On Return ! ! A an upper triangular matrix and the multipliers ! which were used to obtain it. ! The factorization can be written A = L*U where ! L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! IPVT INTEGER(N) ! an INTEGER vector of pivot indices. ! ! RCOND DOUBLE PRECISION ! an estimate of the reciprocal condition of A . ! For the system A*X = B , relative perturbations ! in A and B of size EPSILON may cause ! relative perturbations in X of size EPSILON/RCOND . ! If RCOND is so small that the logical expression ! 1.0 + RCOND .EQ. 1.0 ! is true, then A may be singular to working ! precision. In particular, RCOND is zero if ! exact singularity is detected or the estimate ! underflows. ! ! Z DOUBLE PRECISION(N) ! a work vector whose contents are usually unimportant. ! If A is close to a singular matrix, then Z is ! an approximate null vector in the sense that ! NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . ! !***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. ! Stewart, LINPACK Users' Guide, SIAM, 1979. !***ROUTINES CALLED DASUM, DAXPY, DDOT, DGEFA, DSCAL !***REVISION HISTORY (YYMMDD) ! 780814 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DGECO IMPLICIT NONE INTEGER LDA,N,IPVT(*) DOUBLE PRECISION A(LDA,*),Z(*) DOUBLE PRECISION RCOND ! DOUBLE PRECISION DDOT,EK,T,WK,WKM DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L ! ! COMPUTE 1-NORM OF A ! !***FIRST EXECUTABLE STATEMENT DGECO ANORM = 0.0D0 DO 10 J = 1, N ANORM = MAX(ANORM,DASUM(N,A(1,J),1)) 10 CONTINUE ! ! FACTOR ! CALL DGEFA(A,LDA,N,IPVT,INFO) ! ! RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . ! ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . ! TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE ! CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE ! TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID ! OVERFLOW. ! ! SOLVE TRANS(U)*W = E ! EK = 1.0D0 DO 20 J = 1, N Z(J) = 0.0D0 20 CONTINUE DO 100 K = 1, N IF (Z(K) .NE. 0.0D0) EK = SIGN(EK,-Z(K)) IF (ABS(EK-Z(K)) .LE. ABS(A(K,K))) GO TO 30 S = ABS(A(K,K))/ABS(EK-Z(K)) CALL DSCAL(N,S,Z,1) EK = S*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = ABS(WK) SM = ABS(WKM) IF (A(K,K) .EQ. 0.0D0) GO TO 40 WK = WK/A(K,K) WKM = WKM/A(K,K) GO TO 50 40 CONTINUE WK = 1.0D0 WKM = 1.0D0 50 CONTINUE KP1 = K + 1 IF (KP1 .GT. N) GO TO 90 DO 60 J = KP1, N SM = SM + ABS(Z(J)+WKM*A(K,J)) Z(J) = Z(J) + WK*A(K,J) S = S + ABS(Z(J)) 60 CONTINUE IF (S .GE. SM) GO TO 80 T = WKM - WK WK = WKM DO 70 J = KP1, N Z(J) = Z(J) + T*A(K,J) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) ! ! SOLVE TRANS(L)*Y = W ! DO 120 KB = 1, N K = N + 1 - KB IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0D0) GO TO 110 S = 1.0D0/ABS(Z(K)) CALL DSCAL(N,S,Z,1) 110 CONTINUE L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) ! YNORM = 1.0D0 ! ! SOLVE L*V = Y ! DO 140 K = 1, N L = IPVT(K) T = Z(L) Z(L) = Z(K) Z(K) = T IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) IF (ABS(Z(K)) .LE. 1.0D0) GO TO 130 S = 1.0D0/ABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM ! ! SOLVE U*Z = V ! DO 160 KB = 1, N K = N + 1 - KB IF (ABS(Z(K)) .LE. ABS(A(K,K))) GO TO 150 S = ABS(A(K,K))/ABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 T = -Z(K) CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE ! MAKE ZNORM = 1.0 S = 1.0D0/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM ! IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 RETURN END !*************************************************************************** !*************************************************************************** !DECK DGEDI SUBROUTINE DGEDI (A, LDA, N, IPVT, DET, WORK, JOB) !***BEGIN PROLOGUE DGEDI !***PURPOSE Compute the determinant and inverse of a matrix using the ! factors computed by DGECO or DGEFA. !***LIBRARY SLATEC (LINPACK) !***CATEGORY D3A1, D2A1 !***TYPE DOUBLE PRECISION (SGEDI-S, DGEDI-D, CGEDI-C) !***KEYWORDS DETERMINANT, INVERSE, LINEAR ALGEBRA, LINPACK, MATRIX !***AUTHOR Moler, C. B., (U. of New Mexico) !***DESCRIPTION ! ! DGEDI computes the determinant and inverse of a matrix ! using the factors computed by DGECO or DGEFA. ! ! On Entry ! ! A DOUBLE PRECISION(LDA, N) ! the output from DGECO or DGEFA. ! ! LDA INTEGER ! the leading dimension of the array A . ! ! N INTEGER ! the order of the matrix A . ! ! IPVT INTEGER(N) ! the pivot vector from DGECO or DGEFA. ! ! WORK DOUBLE PRECISION(N) ! work vector. Contents destroyed. ! ! JOB INTEGER ! = 11 both determinant and inverse. ! = 01 inverse only. ! = 10 determinant only. ! ! On Return ! ! A inverse of original matrix if requested. ! Otherwise unchanged. ! ! DET DOUBLE PRECISION(2) ! determinant of original matrix if requested. ! Otherwise not referenced. ! Determinant = DET(1) * 10.0**DET(2) ! with 1.0 .LE. ABS(DET(1)) .LT. 10.0 ! or DET(1) .EQ. 0.0 . ! ! Error Condition ! ! A division by zero will occur if the input factor contains ! a zero on the diagonal and the inverse is requested. ! It will not occur if the subroutines are called correctly ! and if DGECO has set RCOND .GT. 0.0 or DGEFA has set ! INFO .EQ. 0 . ! !***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. ! Stewart, LINPACK Users' Guide, SIAM, 1979. !***ROUTINES CALLED DAXPY, DSCAL, DSWAP !***REVISION HISTORY (YYMMDD) ! 780814 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DGEDI IMPLICIT NONE INTEGER LDA,N,IPVT(*),JOB DOUBLE PRECISION A(LDA,*),DET(2),WORK(*) ! DOUBLE PRECISION T DOUBLE PRECISION TEN INTEGER I,J,K,KB,KP1,L,NM1 !***FIRST EXECUTABLE STATEMENT DGEDI ! ! COMPUTE DETERMINANT ! IF (JOB/10 .EQ. 0) GO TO 70 DET(1) = 1.0D0 DET(2) = 0.0D0 TEN = 10.0D0 DO 50 I = 1, N IF (IPVT(I) .NE. I) DET(1) = -DET(1) DET(1) = A(I,I)*DET(1) IF (DET(1) .EQ. 0.0D0) GO TO 60 10 IF (ABS(DET(1)) .GE. 1.0D0) GO TO 20 DET(1) = TEN*DET(1) DET(2) = DET(2) - 1.0D0 GO TO 10 20 CONTINUE 30 IF (ABS(DET(1)) .LT. TEN) GO TO 40 DET(1) = DET(1)/TEN DET(2) = DET(2) + 1.0D0 GO TO 30 40 CONTINUE 50 CONTINUE 60 CONTINUE 70 CONTINUE ! ! COMPUTE INVERSE(U) ! IF (MOD(JOB,10) .EQ. 0) GO TO 150 DO 100 K = 1, N A(K,K) = 1.0D0/A(K,K) T = -A(K,K) CALL DSCAL(K-1,T,A(1,K),1) KP1 = K + 1 IF (N .LT. KP1) GO TO 90 DO 80 J = KP1, N T = A(K,J) A(K,J) = 0.0D0 CALL DAXPY(K,T,A(1,K),1,A(1,J),1) 80 CONTINUE 90 CONTINUE 100 CONTINUE ! ! FORM INVERSE(U)*INVERSE(L) ! NM1 = N - 1 IF (NM1 .LT. 1) GO TO 140 DO 130 KB = 1, NM1 K = N - KB KP1 = K + 1 DO 110 I = KP1, N WORK(I) = A(I,K) A(I,K) = 0.0D0 110 CONTINUE DO 120 J = KP1, N T = WORK(J) CALL DAXPY(N,T,A(1,J),1,A(1,K),1) 120 CONTINUE L = IPVT(K) IF (L .NE. K) CALL DSWAP(N,A(1,K),1,A(1,L),1) 130 CONTINUE 140 CONTINUE 150 CONTINUE RETURN END !*************************************************************************** !*************************************************************************** !DECK DGEFA SUBROUTINE DGEFA (A, LDA, N, IPVT, INFO) !***BEGIN PROLOGUE DGEFA !***PURPOSE Factor a matrix using Gaussian elimination. !***LIBRARY SLATEC (LINPACK) !***CATEGORY D2A1 !***TYPE DOUBLE PRECISION (SGEFA-S, DGEFA-D, CGEFA-C) !***KEYWORDS GENERAL MATRIX, LINEAR ALGEBRA, LINPACK, ! MATRIX FACTORIZATION !***AUTHOR Moler, C. B., (U. of New Mexico) !***DESCRIPTION ! ! DGEFA factors a double precision matrix by Gaussian elimination. ! ! DGEFA is usually called by DGECO, but it can be called ! directly with a saving in time if RCOND is not needed. ! (Time for DGECO) = (1 + 9/N)*(Time for DGEFA) . ! ! On Entry ! ! A DOUBLE PRECISION(LDA, N) ! the matrix to be factored. ! ! LDA INTEGER ! the leading dimension of the array A . ! ! N INTEGER ! the order of the matrix A . ! ! On Return ! ! A an upper triangular matrix and the multipliers ! which were used to obtain it. ! The factorization can be written A = L*U where ! L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! IPVT INTEGER(N) ! an integer vector of pivot indices. ! ! INFO INTEGER ! = 0 normal value. ! = K if U(K,K) .EQ. 0.0 . This is not an error ! condition for this subroutine, but it does ! indicate that DGESL or DGEDI will divide by zero ! if called. Use RCOND in DGECO for a reliable ! indication of singularity. ! !***REFERENCES J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W. ! Stewart, LINPACK Users' Guide, SIAM, 1979. !***ROUTINES CALLED DAXPY, DSCAL, IDAMAX !***REVISION HISTORY (YYMMDD) ! 780814 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900326 Removed duplicate information from DESCRIPTION section. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DGEFA IMPLICIT NONE INTEGER LDA,N,IPVT(*),INFO DOUBLE PRECISION A(LDA,*) ! DOUBLE PRECISION T INTEGER IDAMAX,J,K,KP1,L,NM1 ! ! GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING ! !***FIRST EXECUTABLE STATEMENT DGEFA INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 ! ! FIND L = PIVOT INDEX ! L = IDAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L ! ! ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED ! IF (A(L,K) .EQ. 0.0D0) GO TO 40 ! ! INTERCHANGE IF NECESSARY ! IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE ! ! COMPUTE MULTIPLIERS ! T = -1.0D0/A(K,K) CALL DSCAL(N-K,T,A(K+1,K),1) ! ! ROW ELIMINATION WITH COLUMN INDEXING ! DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. 0.0D0) INFO = N RETURN END !*************************************************************************** !*************************************************************************** subroutine dgesl(a,lda,n,ipvt,b,job) integer lda,n,ipvt(1),job double precision a(lda,1),b(1) ! ! dgesl solves the double precision system ! a * x = b or trans(a) * x = b ! using the factors computed by dgeco or dgefa. ! ! on entry ! ! a double precision(lda, n) ! the output from dgeco or dgefa. ! ! lda integer ! the leading dimension of the array a . ! ! n integer ! the order of the matrix a . ! ! ipvt integer(n) ! the pivot vector from dgeco or dgefa. ! ! b double precision(n) ! the right hand side vector. ! ! job integer ! = 0 to solve a*x = b , ! = nonzero to solve trans(a)*x = b where ! trans(a) is the transpose. ! ! on return ! ! b the solution vector x . ! ! error condition ! ! a division by zero will occur if the input factor contains a ! zero on the diagonal. technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of lda . it will not occur if the subroutines are ! called correctly and if dgeco has set rcond .gt. 0.0 ! or dgefa has set info .eq. 0 . ! ! to compute inverse(a) * c where c is a matrix ! with p columns ! call dgeco(a,lda,n,ipvt,rcond,z) ! if (rcond is too small) go to ... ! do 10 j = 1, p ! call dgesl(a,lda,n,ipvt,c(1,j),0) ! 10 continue ! ! linpack. this version dated 08/14/78 . ! cleve moler, university of new mexico, argonne national lab. ! ! subroutines and functions ! ! blas daxpy,ddot ! ! internal variables ! double precision ddot,t integer k,kb,l,nm1 ! nm1 = n - 1 if (job .ne. 0) go to 50 ! ! job = 0 , solve a * x = b ! first solve l*y = b ! if (nm1 .lt. 1) go to 30 do 20 k = 1, nm1 l = ipvt(k) t = b(l) if (l .eq. k) go to 10 b(l) = b(k) b(k) = t 10 continue call daxpy(n-k,t,a(k+1,k),1,b(k+1),1) 20 continue 30 continue ! ! now solve u*x = y ! do 40 kb = 1, n k = n + 1 - kb b(k) = b(k)/a(k,k) t = -b(k) call daxpy(k-1,t,a(1,k),1,b(1),1) 40 continue go to 100 50 continue ! ! job = nonzero, solve trans(a) * x = b ! first solve trans(u)*y = b ! do 60 k = 1, n t = ddot(k-1,a(1,k),1,b(1),1) b(k) = (b(k) - t)/a(k,k) 60 continue ! ! now solve trans(l)*x = y ! if (nm1 .lt. 1) go to 90 do 80 kb = 1, nm1 k = n - kb b(k) = b(k) + ddot(n-k,a(k+1,k),1,b(k+1),1) l = ipvt(k) if (l .eq. k) go to 70 t = b(l) b(l) = b(k) b(k) = t 70 continue 80 continue 90 continue 100 continue return end !*************************************************************************** !*************************************************************************** !DECK DSCAL SUBROUTINE DSCAL (N, DA, DX, INCX) !***BEGIN PROLOGUE DSCAL !***PURPOSE Multiply a vector by a constant. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A6 !***TYPE DOUBLE PRECISION (SSCAL-S, DSCAL-D, CSCAL-C) !***KEYWORDS BLAS, LINEAR ALGEBRA, SCALE, VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DA double precision scale factor ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! ! --Output-- ! DX double precision result (unchanged if N.LE.0) ! ! Replace double precision DX by double precision DA*DX. ! For I = 0 to N-1, replace DX(IX+I*INCX) with DA * DX(IX+I*INCX), ! where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900821 Modified to correct problem with a negative increment. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DSCAL IMPLICIT NONE DOUBLE PRECISION DA, DX(*) INTEGER I, INCX, IX, M, MP1, N !***FIRST EXECUTABLE STATEMENT DSCAL IF (N .LE. 0) RETURN IF (INCX .EQ. 1) GOTO 20 ! ! Code for increment not equal to 1. ! IX = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 DO 10 I = 1,N DX(IX) = DA*DX(IX) IX = IX + INCX 10 CONTINUE RETURN ! ! Code for increment equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 5. ! 20 M = MOD(N,5) IF (M .EQ. 0) GOTO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF (N .LT. 5) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I+1) = DA*DX(I+1) DX(I+2) = DA*DX(I+2) DX(I+3) = DA*DX(I+3) DX(I+4) = DA*DX(I+4) 50 CONTINUE RETURN END !*************************************************************************** !*************************************************************************** !DECK DSWAP SUBROUTINE DSWAP (N, DX, INCX, DY, INCY) !***BEGIN PROLOGUE DSWAP !***PURPOSE Interchange two vectors. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A5 !***TYPE DOUBLE PRECISION (SSWAP-S, DSWAP-D, CSWAP-C, ISWAP-I) !***KEYWORDS BLAS, INTERCHANGE, LINEAR ALGEBRA, VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! DY double precision vector with N elements ! INCY storage spacing between elements of DY ! ! --Output-- ! DX input vector DY (unchanged if N .LE. 0) ! DY input vector DX (unchanged if N .LE. 0) ! ! Interchange double precision DX and double precision DY. ! For I = 0 to N-1, interchange DX(LX+I*INCX) and DY(LY+I*INCY), ! where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is ! defined in a similar way using INCY. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890831 Modified array declarations. (WRB) ! 890831 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 920310 Corrected definition of LX in DESCRIPTION. (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE DSWAP IMPLICIT NONE INTEGER N,INCX,INCY,IX,IY,I,M,MP1,NS DOUBLE PRECISION DX(*), DY(*), DTEMP1, DTEMP2, DTEMP3 !***FIRST EXECUTABLE STATEMENT DSWAP IF (N .LE. 0) RETURN IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60 ! ! Code for unequal or nonpositive increments. ! 5 IX = 1 IY = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 IF (INCY .LT. 0) IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP1 = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP1 IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN ! ! Code for both increments equal to 1. ! ! Clean-up loop so remaining vector length is a multiple of 3. ! 20 M = MOD(N,3) IF (M .EQ. 0) GO TO 40 DO 30 I = 1,M DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 30 CONTINUE IF (N .LT. 3) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP1 = DX(I) DTEMP2 = DX(I+1) DTEMP3 = DX(I+2) DX(I) = DY(I) DX(I+1) = DY(I+1) DX(I+2) = DY(I+2) DY(I) = DTEMP1 DY(I+1) = DTEMP2 DY(I+2) = DTEMP3 50 CONTINUE RETURN ! ! Code for equal, positive, non-unit increments. ! 60 NS = N*INCX DO 70 I = 1,NS,INCX DTEMP1 = DX(I) DX(I) = DY(I) DY(I) = DTEMP1 70 CONTINUE RETURN END !*************************************************************************** !*************************************************************************** !DECK IDAMAX INTEGER FUNCTION IDAMAX (N, DX, INCX) !***BEGIN PROLOGUE IDAMAX !***PURPOSE Find the smallest index of that component of a vector ! having the maximum magnitude. !***LIBRARY SLATEC (BLAS) !***CATEGORY D1A2 !***TYPE DOUBLE PRECISION (ISAMAX-S, IDAMAX-D, ICAMAX-C) !***KEYWORDS BLAS, LINEAR ALGEBRA, MAXIMUM COMPONENT, VECTOR !***AUTHOR Lawson, C. L., (JPL) ! Hanson, R. J., (SNLA) ! Kincaid, D. R., (U. of Texas) ! Krogh, F. T., (JPL) !***DESCRIPTION ! ! B L A S Subprogram ! Description of Parameters ! ! --Input-- ! N number of elements in input vector(s) ! DX double precision vector with N elements ! INCX storage spacing between elements of DX ! ! --Output-- ! IDAMAX smallest index (zero if N .LE. 0) ! ! Find smallest index of maximum magnitude of double precision DX. ! IDAMAX = first I, I = 1 to N, to maximize ABS(DX(IX+(I-1)*INCX)), ! where IX = 1 if INCX .GE. 0, else IX = 1+(1-N)*INCX. ! !***REFERENCES C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. ! Krogh, Basic linear algebra subprograms for Fortran ! usage, Algorithm No. 539, Transactions on Mathematical ! Software 5, 3 (September 1979), pp. 308-323. !***ROUTINES CALLED (NONE) !***REVISION HISTORY (YYMMDD) ! 791001 DATE WRITTEN ! 890531 Changed all specific intrinsics to generic. (WRB) ! 890531 REVISION DATE from Version 3.2 ! 891214 Prologue converted to Version 4.0 format. (BAB) ! 900821 Modified to correct problem with a negative increment. ! (WRB) ! 920501 Reformatted the REFERENCES section. (WRB) !***END PROLOGUE IDAMAX IMPLICIT NONE DOUBLE PRECISION DX(*), DMAX, XMAG INTEGER I, INCX, IX, N !***FIRST EXECUTABLE STATEMENT IDAMAX IDAMAX = 0 IF (N .LE. 0) RETURN IDAMAX = 1 IF (N .EQ. 1) RETURN ! IF (INCX .EQ. 1) GOTO 20 ! ! Code for increments not equal to 1. ! IX = 1 IF (INCX .LT. 0) IX = (-N+1)*INCX + 1 DMAX = ABS(DX(IX)) IX = IX + INCX DO 10 I = 2,N XMAG = ABS(DX(IX)) IF (XMAG .GT. DMAX) THEN IDAMAX = I DMAX = XMAG ENDIF IX = IX + INCX 10 CONTINUE RETURN ! ! Code for increments equal to 1. ! 20 DMAX = ABS(DX(1)) DO 30 I = 2,N XMAG = ABS(DX(I)) IF (XMAG .GT. DMAX) THEN IDAMAX = I DMAX = XMAG ENDIF 30 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MATIDENT(N,IDENTITY) IMPLICIT NONE INTEGER I,J,N DOUBLE PRECISION IDENTITY(N,N) ! PRODUCE AN IDENTITY MATRIX OF SIZE NxN DO 20 I=1,N DO 10 J=1,N IF (I.EQ.J) THEN IDENTITY(I,J) = 1.0D0 ELSE IDENTITY(I,J) = 0.0D0 END IF 10 CONTINUE 20 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MATDIAG(N,VEC,DIAG) IMPLICIT NONE INTEGER I,J,N DOUBLE PRECISION VEC(N),DIAG(N,N) ! PRODUCE A DIAGONAL MATRIX OF SIZE NxN DO 20 I=1,N DO 10 J=1,N IF (I.EQ.J) THEN DIAG(I,J) = VEC(I) ELSE DIAG(I,J) = 0.0D0 END IF 10 CONTINUE 20 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MM_DG(N,DVEC,G,PROD) IMPLICIT NONE INTEGER I,J,N DOUBLE PRECISION DVEC(N),G(N,N),PROD(N,N) ! MULTIPLY GENERAL MATRIX G BY DIAGONAL MATRIX D ON THE LEFT ! TO PRODUCE MATRIX PROD DO 20 I=1,N DO 10 J=1,N PROD(I,J) = DVEC(I)*G(I,J) 10 CONTINUE 20 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MM_D1GD2(N,D1VEC,G,D2VEC,PROD) IMPLICIT NONE INTEGER I,J,N DOUBLE PRECISION D1VEC(N),D2VEC(N),G(N,N),PROD(N,N) ! MULTIPLY GENERAL MATRIX G BY DIAGONAL MATRIX D1 ON THE LEFT ! AND DIAGONAL MATRIX D2 ON THE RIGHT TO PRODUCE MATRIX PROD DO 20 I=1,N DO 10 J=1,N PROD(I,J) = D1VEC(I)*G(I,J)*D2VEC(J) 10 CONTINUE 20 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MM_GD(N,G,DVEC,PROD) IMPLICIT NONE INTEGER I,J,N DOUBLE PRECISION DVEC(N),G(N,N),PROD(N,N) ! MULTIPLY GENERAL MATRIX G BY DIAGONAL MATRIX D ON THE RIGHT ! TO PRODUCE MATRIX PROD DO 20 J=1,N DO 10 I=1,N PROD(I,J) = G(I,J)*DVEC(J) 10 CONTINUE 20 CONTINUE RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MM_IG1G2(N,M,G1,G2,PROD,INV_IO) !INPUT : N,M,G1,G2,INV_IO !OUTPUT: PROD !THIS PROGRAM COMPUTES THE PRODUCT OF THE INVERSE OF A GENERAL ! (NONSINGULAR) NxN MATRIX G1 AND A GENERAL NxM MATRIX G2 !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY MM_IG1G2********************************** ! MATMUL !************************************************************************ !EXTERNAL SUBPROGRAMS USED BY MM_IG1G2*********************************** ! DGEFA,DGECO,DGESL !************************************************************************ IMPLICIT NONE !INPUT VARIABLES INTEGER :: & N,M,INV_IO DOUBLE PRECISION, DIMENSION(N,N) :: & G1 DOUBLE PRECISION, DIMENSION(N,M) :: & G2 !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N,M) :: & PROD !INTERNAL VARIABLES INTEGER :: & I,J,INFO,JOB INTEGER, DIMENSION(N) :: & IPVT DOUBLE PRECISION :: & RCOND DOUBLE PRECISION, DIMENSION(N) :: & Z DOUBLE PRECISION, DIMENSION(N,N) :: & G1TEMP DOUBLE PRECISION, DIMENSION(N,M) :: & G2TEMP !COMPUTE ONLY THE L*U FACTORIZATION OF MATRIX G1 FOR SUBROUTINE DGESL ! IF NOT DOING DIAGNOSTIC TESTING (I.E. NORMAL OPERATION) IF (INV_IO == 0) THEN G1TEMP = G1 CALL DGEFA(G1TEMP,N,N,IPVT,INFO) !COMPUTE THE L*U FACTORIZATION OF MATRIX G1 FOR SUBROUTINE !DGESL AND ITS RECIPROCAL CONDITION NUMBER FOR DIAGNOSTIC !TESTING PURPOSES ELSE G1TEMP = G1 CALL DGECO(G1TEMP,N,N,IPVT,RCOND,Z) PRINT* PRINT*,'THE RECIPROCAL CONDITION NUMBER OF THE MATRIX IS:' PRINT*, RCOND END IF !COMPUTE THE PRODUCT OF THE INVERSE OF MATRIX G1 WITH MATRIX G2 USING ! OUTPUT OF SUBROUTINE DGEFA OR DGECO G2TEMP = G2 DO J=1,M JOB = 0 CALL DGESL(G1TEMP,N,N,IPVT,G2TEMP(1,J),JOB) END DO PROD = G2TEMP IF (INV_IO == 1) THEN PRINT* PRINT*,'THE MATRIX PRODUCT LOOKS LIKE:' WRITE(*,10) ((PROD(I,J),J=1,M),I=1,N) 10 FORMAT(8(1X,F14.8)) END IF RETURN END !************************************************************************ !************************************************************************ SUBROUTINE MV_DV(N,DVEC,V,PROD) IMPLICIT NONE INTEGER I,N DOUBLE PRECISION DVEC(N),V(N),PROD(N) ! MULTIPLY VECTOR V BY DIAGONAL MATRIX D ! TO PRODUCE MATRIX PROD DO 10 I=1,N PROD(I) = DVEC(I)*V(I) 10 CONTINUE RETURN END !********************************************************************* !********************************************************************* SUBROUTINE PLEG(ANGLEO,GMU,M,NEXP,NZEN,YPLEG) !----------------------------------------------------------------------C ! Computation of renormalized associated Legendre polynomials C ! of the form: C ! C ! m 1/2 m C ! Y (u) = [(l-m)!/(l+m)!] P (u) C ! l l C ! C ! where, C ! m C ! P (u) = Associated Legendre polynomial of C ! l order l and degree m C ! C ! u = Cosine of zenith angle C ! C ! Reference: C ! C ! Dave, J. V.,and Armstrong, B. H., 1970: Computations C ! of High-order Associated Legendre Polynomials, C ! JQSRT, 10, 557-562, 1970 C !----------------------------------------------------------------------C ! I N P U T V A R I A B L E S : C !----------------------------------------------------------------------C ! ANGLEO : Cosine of solar zenith angle C ! GMU : Gaussian quadrature points C ! M : Index for degree of Legendre polynomial C ! NEXP : Tot. no. of polynomial expansion terms C ! NZEN : Tot. no. of quadrature points C !----------------------------------------------------------------------C ! O U T P U T V A R I A B L E S : C !----------------------------------------------------------------------C ! YPLEG : Renormalized associated Legendre polynomials C !----------------------------------------------------------------------C IMPLICIT NONE INTEGER NEXP,NZEN,M,I,IL,L DOUBLE PRECISION CM1,DM1,CM2,DM2 DOUBLE PRECISION ANGLEO,GMU(NZEN),LM1,LM2,MU,YPLEG(NEXP,NZEN+2) SAVE CM1, DM1 ! IF (M .EQ. 0) THEN CM1 = 1. DM1 = 1. ELSE CM2 = -SQRT((2.*M-1.)/(2.*M))*CM1 DM2 = -SQRT((2.*M+1.)/(2.*M))*DM1 END IF DO 20,I=1,NZEN+2 IF (I .LE. NZEN) MU = GMU(I) IF (I .EQ. NZEN+1) MU = ANGLEO IF (I .EQ. NZEN+2) MU = 0. !----------------------------------------------------------------------C ! Compute initial values C !----------------------------------------------------------------------C IF (M .EQ. 0) THEN YPLEG(1,I) = 1. YPLEG(2,I) = MU ELSE YPLEG(M+1,I) = CM2*(1-MU**2)**(M/2.) YPLEG(M+2,I) = DM2*MU*(1-MU**2)**(M/2.) END IF DO 10,IL=1,NEXP-2 L = IL-1 IF (IL .LT. M+1) THEN YPLEG(IL,I) = 0. GO TO 10 ELSE !----------------------------------------------------------------------C ! Use varying-degree recurrence formula to compute C ! successive values C !----------------------------------------------------------------------C LM2 = (L-M+2)*(L+M+2) LM1 = (L+M+1)*(L-M+1) YPLEG(IL+2,I) = ((2*L+3)*MU*YPLEG(IL+1,I)-SQRT(LM1)* & YPLEG(IL,I))/SQRT(LM2) END IF 10 CONTINUE 20 CONTINUE IF (M .GT. 0) THEN CM1 = CM2 DM1 = DM2 END IF RETURN END !************************************************************************** !************************************************************************** SUBROUTINE PLEG2(ANGLEO,GMU,M,NEXP,NZEN,YPLEG) !----------------------------------------------------------------------C ! Computation of renormalized associated Legendre polynomials C ! of the form: C ! C ! m 1/2 m C ! Y (u) = [(l-m)!/(l+m)!] P (u) C ! l l C ! C ! where, C ! m C ! P (u) = Associated Legendre polynomial of C ! l order l and degree m C ! C ! u = Cosine of zenith angle C ! C ! Reference: C ! C ! Dave, J. V.,and Armstrong, B. H., 1970: Computations C ! of High-order Associated Legendre Polynomials, C ! JQSRT, 10, 557-562, 1970 C !----------------------------------------------------------------------C ! I N P U T V A R I A B L E S : C !----------------------------------------------------------------------C ! ANGLEO : Cosine of solar zenith angle C ! GMU : Gaussian quadrature points C ! M : Index for degree of Legendre polynomial C ! NEXP : Tot. no. of polynomial expansion terms C ! NZEN : Tot. no. of quadrature points C !----------------------------------------------------------------------C ! O U T P U T V A R I A B L E S : C !----------------------------------------------------------------------C ! YPLEG : Renormalized associated Legendre polynomials C !----------------------------------------------------------------------C IMPLICIT NONE !INPUT VARIABLES INTEGER, INTENT(IN) :: & M,NEXP,NZEN DOUBLE PRECISION, INTENT(IN) :: & ANGLEO DOUBLE PRECISION, DIMENSION(NZEN), INTENT(IN) :: & GMU !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(NEXP,NZEN+2), INTENT(OUT) :: & YPLEG !INTERNAL VARIABLES INTEGER :: & I,IL,L DOUBLE PRECISION :: & CM2,DM2,LM1,LM2,MU DOUBLE PRECISION, SAVE :: & CM1,DM1 !START PROGRAM YPLEG = 0.0D0 !PRINT* !PRINT*,'INSIDE PLEG2: COMING IN' !PRINT*,'THE YPLEGP(L,MU,DEGREE) INSIDE PLEG2 IS FOR DEGREE = :', M !WRITE(*,40) ((YPLEG(L,I),I=1,NZEN+2),L=1,NEXP) IF (M == 0) THEN CM1 = 1.0D0 DM1 = 1.0D0 ELSE CM2 = -DSQRT(((2.0D0*DBLE(M)) - 1.0D0)/(2.0D0*DBLE(M)))*CM1 DM2 = -DSQRT(((2.0D0*DBLE(M)) + 1.0D0)/(2.0D0*DBLE(M)))*DM1 END IF DO I=1,NZEN+2 IF (I <= NZEN) MU = GMU(I) IF (I == NZEN+1) MU = ANGLEO IF (I == NZEN+2) MU = 0.0D0 !----------------------------------------------------------------------C ! Compute initial values C !----------------------------------------------------------------------C IF (M == 0) THEN YPLEG(1,I) = 1.0D0 YPLEG(2,I) = MU ELSE YPLEG(M+1,I) = CM2*(1.0D0 - MU**2)**(DBLE(M)/2.0D0) YPLEG(M+2,I) = DM2*MU*(1.0D0 - MU**2)**(DBLE(M)/2.0D0) END IF DO IL=1,NEXP-2 L = IL-1 IF (IL < M+1) THEN YPLEG(IL,I) = 0.0D0 ELSE !----------------------------------------------------------------------C ! Use varying-degree recurrence formula to compute C ! successive values C !----------------------------------------------------------------------C LM2 = DBLE( (L-M+2)*(L+M+2) ) LM1 = DBLE( (L+M+1)*(L-M+1) ) YPLEG(IL+2,I) = ( (2.0D0*DBLE(L) + 3.0D0)*MU*YPLEG(IL+1,I) & - DSQRT(LM1)*YPLEG(IL,I) )/DSQRT(LM2) END IF END DO END DO IF (M > 0) THEN CM1 = CM2 DM1 = DM2 END IF !PRINT* !PRINT*,'INSIDE PLEG2: GOING OUT' !PRINT*,'THE YPLEGP(L,MU,DEGREE) INSIDE PLEG2 IS FOR DEGREE = :', M !WRITE(*,40) ((YPLEG(L,I),I=1,NZEN+2),L=1,NEXP) !40 FORMAT(7(1X,F17.13)) END SUBROUTINE PLEG2 !************************************************************************** !************************************************************************** SUBROUTINE SYSSOL(N,M,A,B,X,SYS_IO) !INPUT : N,M,A,B,SYS_IO !OUTPUT: X !THIS PROGRAM SOLVES LINEAR SYSTEM(S) OF EQUATIONS VIA GAUSSIAN ! ELIMINATION !PROGRAMMER: MATT CHRISTI !INTRINSIC SUBPROGRAMS USED BY SYSSOL*********************************** ! MATMUL !*********************************************************************** !EXTERNAL SUBPROGRAMS USED BY SYSSOL************************************ ! DGEFA,DGECO,DGESL !*********************************************************************** IMPLICIT NONE !INPUT VARIABLES INTEGER M,N,SYS_IO DOUBLE PRECISION, DIMENSION(N,M) :: & B DOUBLE PRECISION, DIMENSION(N,N) :: & A !OUTPUT VARIABLES DOUBLE PRECISION, DIMENSION(N,M) :: & X !INTERNAL VARIABLES INTEGER I,J,INFO,JOB INTEGER, DIMENSION(N) :: & IPVT DOUBLE PRECISION RCOND DOUBLE PRECISION DET(2) DOUBLE PRECISION, DIMENSION(N) :: & Z,WORK,PROD DOUBLE PRECISION, DIMENSION(N,N) :: & A_LU !START PROGRAM IF(SYS_IO == 1) THEN PRINT* PRINT*,'ENTERING SYSSOL' END IF !COMPUTE ONLY THE L*U FACTORIZATION OF MATRIX A FOR SUBROUTINE DGESL ! IF NOT DOING DIAGNOSTIC TESTING (I.E. NORMAL OPERATION) IF (SYS_IO == 0) THEN A_LU = A CALL DGEFA(A_LU,N,N,IPVT,INFO) !COMPUTE THE L*U FACTORIZATION OF MATRIX A FOR SUBROUTINE !DGESL AND ITS RECIPROCAL CONDITION NUMBER FOR DIAGNOSTIC !TESTING PURPOSES ELSE A_LU = A CALL DGECO(A_LU,N,N,IPVT,RCOND,Z) PRINT* PRINT*,'THE RECIPROCAL CONDITION NUMBER OF THE MATRIX IS:' PRINT*, RCOND END IF !SOLVE THE LINEAR SYSTEM(S) A*X=B USING OUTPUT OF SUBROUTINE DGEFA OR ! DGECO X = B DO J=1,M JOB = 0 CALL DGESL(A_LU,N,N,IPVT,X(1,J),JOB) IF (SYS_IO == 1) THEN PRINT* PRINT*,'THE SOLUTION X LOOKS LIKE:' WRITE(*,*) (X(I,J),I=1,N) !10 FORMAT(1X,E17.10) PROD = MATMUL(A,X(:,J)) PRINT* PRINT*,'THE PRODUCT OF THE ORIGINAL MATRIX' PRINT*,'AND THE SOLUTION LOOKS LIKE (SHOULD BE = B):' WRITE(*,*) (PROD(I),I=1,N) END IF END DO !END PROGRAM IF (SYS_IO == 1) THEN PRINT* PRINT*,'LEAVING SYSSOL' END IF RETURN END !*********************************************************************** !***********************************************************************