#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'READ.ME' <<'END_OF_FILE' X *************************************************************************** X * All the software contained in this library is protected by copyright. * X * Permission to use, copy, modify, and distribute this software for any * X * purpose without fee is hereby granted, provided that this entire notice * X * is included in all copies of any software which is or includes a copy * X * or modification of this software and in all copies of the supporting * X * documentation for such software. * X *************************************************************************** X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * X * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * X * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * X * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * X * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * X *************************************************************************** X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * X * ABOVE STATEMENT. * X *************************************************************************** X X AUTHORS: X X F. Stenger, B. Keys, M. O'Reilly, K. Parker X University of Utah, Salt Lake City, Utah X X S.-A. Gustafson X Stavanger College, Stavanger, Norway X X X REFERENCE: X X - OPE-IVP-PACK via Sinc indefinite integration and Newton's method X NUMERICAL ALGORITHMS, 20 (1999), PP. 241-268 X X SOFTWARE REVISION DATE: X X MARCH 15, 1999 X X SOFTWARE LANGUAGE: X X FORTRAN 77 X X *************************************************************************** X X XU S I N G T H E P A C K A G E X XPURPOSE XThe package is used to solve numerically an initial-value problem Xfor a system of ordinary differential equations, implementing the Xmethods described in the paper by Stenger et. al., ODE-IVP-PACK Xvia Sinc indefinite integration and Newton's method. X X XCONTENTS XThe package consists of the following items: X This READ.ME file describing how to use the package. See also Section 6 Xin the paper by F. Stenger et. al. X A makefile useful for obtaining the executable code in a Unix system X Files map.h and workspace.h needed to pass the functions associated Xwith the conformal maps and the common block workspace between the fixed Xpart of the package and the file provided by the user, containing the Xproblem to be solved. X File code.f which contains the driver and all the modules that do not need Xbe modified by the user. X 13 further files with extension f, documenting examples in Section 7 Xof the paper which may be run using the package. The name of the files Xand the corresponding problems are given in the table below. XExample 7.1 problem_01.f XExample 7.2 problem_AB.f XExample 7.3 problem_Bul.f XExample 7.4 problem_Fun.f XExample 7.5 problem_Id.f XExample 7.6 problem_N11.f XExample 7.7 problem_RR.f XExample 7.8 problem_Sec.f XExample 7.9 problem_ig.f XExample 7.10 problem_g.f XExample 7.11 problem_sg1.f XExample 7.12 problem_l.f XExample 7.13 problem_wg.f X XSOLVING AN INITIAL-VALUE PROBLEM WITH THE PACKAGE XThe following 5 files are required to produce an executable code: Xmap.h Xworkspace.h Xmakefile Xcode.f Xproblem.f XThe last file defines the example to be solved using the package. XWe will describe below how to write this file X X XSOLVING ON OF THE EXAMPLES IN SECTION 7 OF THE PAPER BY STENGER et. al. XWe discuss first the situation when the user wants to reproduce one of the X13 examples given in Section 7. It will become apparent from what is said Xbelow how to do, if one wants to change one or more of the parameters in Xthese examples. XAssume that the user wants to treat the Problem 7.3 which is contained in Xthe file problem_Bul.f and write the result on the file answer. Then (s)he Xtypes the following commands: X Xcp problem_Bul.f problem.f Xmake Xproblem>answer X XThus the system of ordinary differential equations to be solved is Xspecified on a file which is copied on the file problem.f which Xserves as input for the main program driver. X X In the case of a non-Unix system, the user needs to combine a file Xcontaining source modules, defining the problem to be solved, with the file Xcode.f, then compile, link and run. X X By studying the 13 solved Xexamples the user may get some ideas about how to treat her/his own problem. XThe package may be used either for solving a given problem or for Xverifying that a solution obtained by other means, is correct. X XIn the Examples 7.2, 7.8 and 7.9 not all of the components of the solution Xvector may be expressed as simple functions. In these cases the corresponding Xprograms deliver the observed error only for the simple components. For Xthe others the exact values are replaced by zero and the absolute value Xof the calculated solution appears in the error column of the tables. The Xoutput for these components is not printed in the paper but will be Xdelivered when these programs are run. As well-known many special Xfunctions and families of orthogonal functions are defined as solutions of Xordinary differential equations. See e. g. Examples 7.8 and 7.13 in the paper. XThen this package could be used to tabulate these functions X XSPECIFYING A PROBLEM: XThe user needs to supply the following information in some routines Xcontained in the file problem.f X1. The number of points to use in solving the equation X2. The system of first order differential equations X3. The initial conditions for the system of first order differential equations X4. The Jacobian associated with this system X5. An approximate solution to be used for starting the Newton iterations X6. The conformal map to use X7. The points where the calculated solution should be tabulated X8. The control variable Itest to indicate if the purpose is to verify Xthe correctness of a suggested solution X XIf the system of ordinary differential equations is linear, then the XNewton process converges after one iteration and the starting Xapproximation may be chosen arbitrarily. There is no guarantee that Xthis problem will converge for a nonlinear problem and an error-exit Xhas been provided if an answer is not found after the maximum number of Xiterations indicated by the user X XPROGRAMS NEEDED TO DEFINE A PARTICULAR PROBLEM XThey form the file problem.f and are contained in all of the files Xin the group above. Note that in the sample problems above the X4 subprograms Phi, OneOverPrime, Rho and Psi are defined via Xmap.h XControlVars Xf XJ_ij XPhi XOneOverPhiPrime XRho XPsi XSetXs XSolution XStartAppr XInitCond XAll of these 11 subprogram must be written by the user. X X X XAN EXAMPLE XHere we show how to provide this information for the Example 7.3 Xin Section 7. The corresponding code is in the file problem_Bul.f XIt contains the subprograms: X XSubroutine ControlVars which sets the following variables: X NumEqs, number of equations in the system X SincPts, number of sinc-points X Nxs, number of points in table X MaxIter, maximum number of iterations X Itest, Itest=0, no exact solution known, Itest=1, exact solution known X and is to be verified numerically X XSubroutine F, returns the derivatives of the solution dy(i)/dt X XSubroutine J_ij, returns the derivatives in the Jacobian, dYDOT(i)/dY(k) X XFour functions defining the conformal mapping, namely XPHi, OneOverPhiPrime, Rho and Psi X XFunction Solution which returns the exact solution, if known. Otherwise the Xuser can define its output to be zero X XSubroutine InitCond returns the initial conditions X XSubroutine StartAppr returns the starting approximate solution for the Newton Xprocess. If the system to be solved is linear, the Newton process converges Xafter one iteration for any starting solution. X X X X X X X X X END_OF_FILE if test 7983 -ne `wc -c <'READ.ME'`; then echo shar: \"'READ.ME'\" unpacked with wrong size! fi # end of 'READ.ME' fi if test -f 'makefile' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'makefile'\" else echo shar: Extracting \"'makefile'\" \(572 characters\) sed "s/^X//" >'makefile' <<'END_OF_FILE' X###################################################################### X# f o r t r a n / m a k e f i l e X# X# This is the makefile for the FORTRAN version of the indefinite X# integration package. X# X###################################################################### X X X################################### X# Executables X################################### X Xproblem : problem.o code.o X f77 -o problem problem.o code.o X X################################### X# Object Files X################################### X Xcode.o : code.f workspace.h Xproblem.o : problem.f workspace.h X X END_OF_FILE if test 572 -ne `wc -c <'makefile'`; then echo shar: \"'makefile'\" unpacked with wrong size! fi # end of 'makefile' fi if test -f 'map.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'map.h'\" else echo shar: Extracting \"'map.h'\" \(1436 characters\) sed "s/^X//" >'map.h' <<'END_OF_FILE' X* m a p . h X X* This module declares all the necessary functions associated with the X* supplied conformal maps. The functions are defined the section X* named map.f in the file code.f X X* The user may define additional conformal maps by writing code X* defining new functions PHi, OneOverPhiPrime, Psi and Rho. In X* this case the instruction INCLUDE map.h is deleted and the X* user's own functions occur instead of the standard ones. X X* M A P L I S T X X DOUBLE PRECISION Phi01, PhiPriRec01, Psi01, Rho01 X EXTERNAL Phi01, PhiPriRec01, Psi01, Rho01 X X DOUBLE PRECISION PhiN11, PhiPriRecN11, PsiN11, RhoN11 X EXTERNAL PhiN11, PhiPriRecN11, PsiN11, RhoN11 X X DOUBLE PRECISION PhiAB, PhiPriRecAB, PsiAB, RhoAB X EXTERNAL PhiAB, PhiPriRecAB, PsiAB, RhoAB X X DOUBLE PRECISION PhiBul, PhiPriRecBul, PsiBul, RhoBul X EXTERNAL PhiBul, PhiPriRecBul, PsiBul, RhoBul X X DOUBLE PRECISION PhiSec, PhiPriRecSec, PsiSec, RhoSec X EXTERNAL PhiSec, PhiPriRecSec, PsiSec, RhoSec X X DOUBLE PRECISION PhiRR, PhiPriRecRR, PsiRR, RhoRR X EXTERNAL PhiRR, PhiPriRecRR, PsiRR, RhoRR X X DOUBLE PRECISION PhiId, PhiPriRecId, PsiId, RhoId X EXTERNAL PhiId, PhiPriRecId, PsiId, RhoId X X DOUBLE PRECISION PhiFun, PhiPriRecFun, PsiFun, RhoFun X EXTERNAL PhiFun, PhiPriRecFun, PsiFun, RhoFun X X X X X X X X X X X X X X X X X X X X X X X X END_OF_FILE if test 1436 -ne `wc -c <'map.h'`; then echo shar: \"'map.h'\" unpacked with wrong size! fi # end of 'map.h' fi if test -f 'workspace.h' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'workspace.h'\" else echo shar: Extracting \"'workspace.h'\" \(2511 characters\) sed "s/^X//" >'workspace.h' <<'END_OF_FILE' X**** w o r k s p a c e . h X X* The file workspace.h should be included in the following modules X* of the user-supplied file problem.f namely SUBROUTINE F, X* SUBROUTINE J_ij and SUBROUTINE StartAppr X X* This header file defines the common block Workspace which is used to X* define a large number of pointers into the working array. X X* Points to the list of the Psi values at the sinc points. X INTEGER PsiPointer X X* Pointer to the values of One-Over-Phi-Prime at the sinc points. X INTEGER OneOverPhiPrimePointer X X* Pointer into the Sigma array. Note that this points to X* Sigma(0), the sigma array actually goes form -N to N. X INTEGER OffSigma X X* Points to a list of pointers (one per equation) for the OldY value X* at the sinc points. X INTEGER OldYPointerPointer X X* Points to the list of Y values at the sinc points. X INTEGER Y0Pointer X X* Points to a temporary matrix used for the linear algebra X* sub-problem AX = B. The size is ( (2N + 1) NumOfEqua )^2 X INTEGER APointer X X* Points to a temporary matrix used for the linear algebra X* sub-problem AX = B. The size is (2N + 1) NumOfEqua X INTEGER BPointer X X* Points to a temporary matrix used to store the pivot values for the X* L-U decomposition of A. Note that the actual decomposition is X* stored in A. The size is (2N + 1) NumOfEqua X INTEGER APivotPointer X X* Points to the list of pointers (one per equation) for the F(t,Y) values X* at the sinc points. X INTEGER FPointerPointer X X* Points to a list of pointers (one per equation) for the BijX values X* at the sinc points. X INTEGER BijPointerPointer X X* Points to a list of pointers (one per equation) for the Y values at X* the sinc points X INTEGER YPointerPointer X X* Points to the list of Phi(X) values at the interpolation points. X INTEGER PhiPointer X X* Pointer to the list of Rho (i.e., Exp( Phi(X) ) ) values at the X* interpolation points. X INTEGER RhoPointer X X* Exp(KH) for each K in [-N,N] X INTEGER OffExpKH X X* Pointer to the list of the Sinc function evaluated at X* equally spaced points. X INTEGER SincPointer X X********************************************************************** X COMMON /Workspace/ X . APointer, APivotPointer, BPointer, BijPointerPointer, X . FPointerPointer, OffExpKH, OffSigma, OldYPointerPointer, X . OneOverPhiPrimePointer, PhiPointer, PsiPointer, RhoPointer, X . SincPointer, YPointerPointer, Y0Pointer X X X X X X X X END_OF_FILE if test 2511 -ne `wc -c <'workspace.h'`; then echo shar: \"'workspace.h'\" unpacked with wrong size! fi # end of 'workspace.h' fi if test -f 'code.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'code.f'\" else echo shar: Extracting \"'code.f'\" \(174456 characters\) sed "s/^X//" >'code.f' <<'END_OF_FILE' X* d r i v e r .f X X* STRUCTURE OF THIS FILE, code.f X* It contains: X* MAIN program driver. X* Several groups of subprograms namely: X* 1. General subprograms written by the authors for solving a X* general system of ordinary differential equations: X* Interpo X* Init1Work X* Init2Work X* Newton X* SincMNRD X* UpDateA X* 2. Functions defining standard conformal mappings and communicated to X* problem.f via map.h: X* Used by file problem_Fun.f for Example 7.4 X* phifun: X* psifun: X* phiprirecfun: X* rhofun: X X* Used by problem_Id.f in Example 7.5 X* phiid: X* psiid: X* phiprirecid: X* rhoid: X X* Used by problem_01.f in Example 7.1 X* phi01: X* psi01: X* phiprirec01: X* rho01: X X* Used by problem_N11.f in Example 7.6 X* phin11: X* psin11: X* phiprirecn11: X* rhon11: X X* Used by: problem_AB.f in Example 7.2 X* problem_ig.f in Example 7.9 X* problem_sg1.f.f in Example 7.11 X* problem_l.f in Example 7.12 X* problem_wg.f in Example 7.13 X* phiab: X* psiab: X* phiprirecab: X* rhoab: X* X* Used by: problem_RR.f in Example 7.7 X* phirr: X* psirr: X* phiprirecrr: X* rhorr: X X* Used by: problem_Bul.f in Example 7.3 X* phibul: X* psibul: X* phiprirecbul: X* rhobul: X* X* Used by: problem_Sec.f in Example 7.8 X* problem_g.f in Example 7.10 X* phisec: X* psisec: X* phiprirecsec: X* rhosec: X* 3. Subprograms from public domain libraries. Complete documentation is X* found in the source code. X* General numerical routines: X* DGEFA from Linpack X* DAXPY from Linpack X* DSCAL from Linpack X* IDAMAX from Linpack X* DGESL from Linpack X* DDOT from Linpack X* DGECO from Linpack X* DGEFA from Linpack X* DASUM form Linpack X* DGAMI from LANL X* DGAMIT from LANL X* DGAMR from LANL X* DGAMMA from LANL X* DGAMLM from LANL X* DLNGAM from LANL X* D1MACH from LANL X* D9LGMC from LANL X* D9GMIT from LANL X* D9LGIC from LANL X* D9LGIT from LANL X* DCSEVL from LANL X* INITDS from LANL X* DLGAMS from LANL X* Routines for handling error messages from SLATEC: X* J4SAVE X* XGETUA X* XERABT X* XERCTL X* XERPRT X* XERROR X* XERRWV X* XERCLR X* XERSAV X* XGETF X* XSETF X*4. Hardware dependent subprograms X* CALJY0 from Argonne National laboratory X* BESJ0 from Argonne National laboratory X* BESY0 from Argonne National laboratory X* FDUMP from SLATEC X* I1MACH from Bell Labs X* 5. Intrinsic functions: X* abs X* aint X* cos X* dabs X* dble X* dexp X* dfloat X* dint X* dlog X* dmax1 X* dmin1 X* dmod X* dsign X* dsin X* dsinh X* dsqrt X* dtanh X* float X* iabs X* len X* log X* log10 X* max0 X* min X* min0 X* mod X* nint X* sin X* sngl X* sqrt X X* This MAIN program calls the routines to solve the IVP and to interpolate. X* If itest=1 we know the exact solution and compare the interpolated X* solution to the exact solution. X* If itest=0, the exact solution is not known and we estimate the error X* to be of the order of magnitude of number Tol, calculated as a local X* variable inside INTEGER FUNCTION Newton and printed by this subprogram X* The calculation of Tol is discussed in Subsection 6.11 of the paper X* "ODE-IVP-PACK via Sinc indefinite integration and Newton's method" X* by F. Stenger, S.-A. Gustafson, B. Keys, M. O'Reilly and K. Parker, X* March 1999 X* ----+----------------------------------------------------------------+ X X PROGRAM driver X* CALLED BY: X* none X* SUBPROGRAMS CALLED: X* ControlVars X* Interpo X* Init1Work X* InitCond X* StartAppr X* Newton X* SetXs X* Solution X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* The function that evaluates the analytic solution. X DOUBLE PRECISION Solution X EXTERNAL Solution X X* The function that computes the newton iterates. X INTEGER Newton X EXTERNAL Newton X X* This section defines constants that are required by the driver X* program. X INTEGER MaxNxs X PARAMETER ( MaxNxs = 500 ) X X INTEGER MaxSize X PARAMETER ( MaxSize = 160000 ) X X INTEGER MaxNumOfEqua X PARAMETER ( MaxNumOfEqua = 8 ) X X* This section allocates memory for all the arrays. X INTEGER N, NumOfEqua, Nxs X INTEGER i, k, Status, TotalAdd, MaxIter, Itest X X DOUBLE PRECISION H X DOUBLE PRECISION TrueValue, Error X DOUBLE PRECISION InitC(MaxNumOfEqua), X(MaxNxs), Y(MaxNxs) X DOUBLE PRECISION Work(MaxSize) X X* ----+----------------------------------------------------------------+ X* This is the executable code to solve the IVP at the Sinc points. X* ----+----------------------------------------------------------------+ X* Input control parameters. X CALL ControlVars( N, NumOfEqua, Nxs, MaxIter, H ,Itest) X X* Next, Initialize the Work array for the chosen conformal mapping. X CALL Init1Work(N, NumOfEqua, H, Work, TotalAdd) X X* Set the initial conditions for the IVP here. X CALL InitCond(NumOfEqua, InitC) X X* Set the initial values for the dependent variables here. X CALL StartAppr( NumOfEqua, N, Work ) X X* CALL the IVP solver. X Status = Newton(H, N, NumOfEqua, InitC, Work, TotalAdd, MaxIter ) X X IF (Status .NE. 0) THEN X PRINT *,' Newton did not converge --> STOP !' X STOP X END IF X X* ----+----------------------------------------------------------------+ X* This is the executable code to interpolate the solution between the X* Sinc points X* ----+----------------------------------------------------------------+ X X* Set k equal to the function you wish to work with. X DO k = 1, NumOfEqua X X PRINT *,' Equation ', k X PRINT *,' ' X IF(Itest.eq.1) THEN X PRINT *,' ','argument, calc. Y_k, exact Y_k',' error' X ELSE X PRINT *,' ',' ',' ',' ',' argument, Y_k-value ' X ENDIF X* Setup a vector of x-values at which you want the solution interpolated. X CALL SetXs( k, Nxs, X ) X X* Call the interpolation routine. X CALL Interpo(H, N, k, Work, Nxs, X, Y, TotalAdd ) X X* Compare analytic solution and computed solution. X IF(Itest .EQ. 1) THEN X DO i = 1, Nxs X TrueValue = Solution(k, X(i)) X Error = DABS(Y(i) - TrueValue) X X write(*,10) X(I), '&',Y(I) ,'&',TrueValue,'&',Error X 10 format(( x, e11.3, x, a1 ), x, e14.8, x, a1 , X + e14.8,x,a1,( x, e11.3 ) ) X X END DO X ELSE X DO i = 1, Nxs X write(*,20) X(I), '&',Y(I) X 20 format(( x, e11.3, x, a1 ), x, e14.8, x, a1) X END DO X END IF X END DO X STOP X END X X X X X SUBROUTINE Interpo(H, N, EquaNum, Work, Nxs, X, Y, TotalAdd) X* CALLED BY: X* driver X* SUBPROGRAMS CALLED: X* Init2Work X* SincMNRD X* i n t e r p o . f X X* This subroutine assumes that newton has already been called and X* that its results have been placed in the work array. X X* ARGUMENT LIST X X* H - DOUBLE PRECISION - This must be the same value as was used on X* the CALL to Newton in order fo the subroutine to correctly X* interpet the values stored in the work array. X X* N - INTEGER - The `size' of the sinc approximation. The total X* number of terms in the approximation is actually 2N + 1. X X* EquaNum - INTEGER - The number of dependent variables to be X* reported. X X* Work - DOUBLE PRECISION array - Containing all intermediate results. X X* Nxs - INTEGER number of independent variable values (e.g., X* x-values) for which the user would like Interpolated values of the X* solution. X X* X(Nxs) - DOUBLE PRECISION array - The independent variable values X* for which the user wants interpolated values of the solution. X X* Y(NumOfEqua, Nxs) - DOUBLE PRECISION array - The interpolated X* values of the dependent variable corresponding to the values X* in X. X X* TotalAdd - INTEGER - Points to the first available free space X* in the work array. X X* CREDITS: X X* This code was written by Michael O'Reilly, Brian Keyes, and X* Kenneth Parker, at the University of Utah, March 20, 1991. It was X* rewritten by Kenneth Parker, Michael O'Reilly, and Brian Keyes, at X* the University of Utah, July 30, 1992. All work was performed X* under the direction of Frank Stenger X X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* This section defines the argument list. X INTEGER N, EquaNum, Nxs, TotalAdd X X DOUBLE PRECISION H X DOUBLE PRECISION Work(*), X(Nxs), Y(Nxs) X X X* This section defines local variables. X INTEGER i, l, ll, YPointer X X DOUBLE PRECISION Check X X* Start of executable code. X l = 2*N+1 X X X* Interpolate the solution at the points X. The assumption is X* that the work array already contains the solution at the sinc points. X X CALL Init2Work( Nxs, X, Work, TotalAdd) X X YPointer = NINT(Work(YPointerPointer + EquaNum)) X DO i = 1, Nxs X Y(i) = (Work(YPointer + 1) + X . Work(RhoPointer + i) * Work(YPointer + 2*N+1)) / X . (1.0D0 + Work(RhoPointer + i)) X CALL SincMNRD(H, -N, N, Work(PhiPointer + i), Work) X X DO ll = 1, 2*N + 1 X Check = Work(YPointer + ll) - X . (Work(YPointer + 1) + X . Work(OffExpKH + ll)*Work(YPointer + 2*N+1)) / X . (1.0D0 + Work(OffExpKH + ll)) X Y(i) = Y(i) + Check * Work(SincPointer + ll) X END DO X END DO X X RETURN X END X X* ----+----------------------------------------------------------------+ X X* m a p . f X X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the X* conformal map from the real line to the real line used to X* approximate functions that decay algebraically at negative X* infinity and exponentially at infinity. X* CALLED BY: problem_Fun.f in Example 7.4 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiFun( z ) X X DOUBLE PRECISION z X X PhiFun = DLOG( DSINH( z + DSQRT( 1.0d0 + z*z ) ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiFun( w ) X X DOUBLE PRECISION w X DOUBLE PRECISION t X X t = dlog( dexp( w ) + dsqrt( 1.0d0 + dexp( 2.0d0 * w ) ) ) X X PsiFun = ( t - 1.0d0/t ) / 2.0d0 X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecFun( z ) X X DOUBLE PRECISION z X DOUBLE PRECISION root X X root = DSQRT( 1.0d0 + z*z ) X X PhiPriRecFun = dtanh( z + root ) * root / ( z + root ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoFun( z ) X X DOUBLE PRECISION z X X RhoFun = dsinh( z + dsqrt( 1.0d0 + z*z ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the identity X* conformal map from the real line to the real line used to X* approximate functions that decay exponentially at negative X* infinity and exponentially at infinity. X* CALLED BY: problem_Id.f in Example 7.5 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiId( z ) X X DOUBLE PRECISION z X X PhiId = z X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiId( w ) X X DOUBLE PRECISION w X X PsiId = w X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecId( z ) X X DOUBLE PRECISION z X X PhiPriRecId = 1.0d0 X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoId( z ) X X DOUBLE PRECISION z X X RhoId = DEXP( z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the conformal X* map from the interval [0,1] to the real line. X* CALLED BY: problem_01.f in Example 7.1 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION Phi01( z ) X X DOUBLE PRECISION z X X Phi01 = DLOG( z / ( 1.0D0 - z ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION Psi01( w ) X X DOUBLE PRECISION w X X Psi01 = 1.0D0 / ( 1.0D0 + DEXP( - w ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRec01( z ) X X DOUBLE PRECISION z X X PhiPriRec01 = z * ( 1.0D0 - z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION Rho01( z ) X X DOUBLE PRECISION z X X Rho01 = z / ( 1.0D0 - z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the conformal X* map from the interval [-1,1] to the real line. X* CALLED BY: problem_N11.f in Example 7.6 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiN11( z ) X X DOUBLE PRECISION z X X PhiN11 = DLOG( ( 1.0D0 + z ) / ( 1.0D0 - z ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiN11( w ) X X DOUBLE PRECISION w X X PsiN11 = DTANH( w / 2.0D0 ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecN11( z ) X X DOUBLE PRECISION z X X PhiPriRecN11 = ( 1.0d0 - z * z ) / 2.0d0 X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoN11( z ) X X DOUBLE PRECISION z X X RhoN11 = ( 1.0d0 + z ) / ( 1.0D0 - z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the conformal X* map from the interval [a,b] to the real line. X* CALLED BY: problem_AB.f in Example 7.2 X* problem_ig.f in Example 7.9 X* problem_sg1.f.f in Example 7.11 X* problem_l.f in Example 7.12 X* problem_wg.f in Example 7.13 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiAB(z,a,b) X X DOUBLE PRECISION a, b, z X X PhiAB = DLOG( ( z - a ) / ( b - z ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiAB(w, a, b) X X DOUBLE PRECISION a, b, w X X PsiAB = (a + b * DEXP(w)) / (1.0D0 + DEXP(w)) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecAB(z, a, b ) X X DOUBLE PRECISION a, b, z X X PhiPriRecAB = (z - a) * ( b - z ) / (b - a) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoAB(z, a, b) X X DOUBLE PRECISION a, b, z X X RhoAB = (z - a) / (b - z) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the X* conformal map from the real line to the real line used to X* approximate functions that decay algebraically at negative X* infinity and algebraically at infinity. X* CALLED BY: problem_RR.f in Example 7.7 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiRR( z ) X X DOUBLE PRECISION z X X PhiRR = DLOG( z + DSQRT( 1.0D0 + z * z ) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiRR( z ) X X DOUBLE PRECISION z X X PsiRR = DSINH( z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecRR( z ) X X DOUBLE PRECISION z X X PhiPriRecRR = DSQRT( 1.0D0 + z * z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoRR( z ) X X DOUBLE PRECISION z X X RhoRR = z + DSQRT( 1.0D0 + z * z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the X* conformal map from the positive real line to the real line used to X* approximate functions that decay algebraically at the origin X* and exponentially at infinity. X* CALLED BY: problem_Bul.f in Example 7.3 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiBul( z ) X X DOUBLE PRECISION z X X PhiBul = DLOG( DSINH(z) ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiBul( w ) X X DOUBLE PRECISION w X X DOUBLE PRECISION a3 X DOUBLE PRECISION a5 X DOUBLE PRECISION a7 X DOUBLE PRECISION a9 X X DOUBLE PRECISION expw X DOUBLE PRECISION expw2 X X a3 = 0.166666666666667D0 X a5 = 0.075000000000000D0 X a7 = 0.044642857142857D0 X a9 = 0.030381944444444 D0 X X expw = DEXP( w ) X expw2 = expw * expw X X IF ( expw .GT. 0.1D0 ) THEN X PsiBul = DLOG( expw + DSQRT( 1.0D0 + expw*expw ) ) X ELSE X PsiBul = ( ( ( ( a9 * expw2 - a7 ) * expw2 + a5 ) * expw2 X . + a3 ) * expw2 + 1.0D0 ) * expw X END IF X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecBul( z ) X X DOUBLE PRECISION z X X PhiPriRecBul = DTANH(z) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoBul( z ) X X DOUBLE PRECISION z X X RhoBul = DSINH(z) X X RETURN X END X* ----+----------------------------------------------------------------+ X* Phi, Psi, 1/Phi' and Rho associated with the X* conformal map from the positive real line to the real line used to X* approximate functions that decay exponentially at the origin X* and exponentially at infinity. X* CALLED BY: problem_Sec.f in Example 7.8 X* problem_g.f in Example 7.10 X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiSec( z ) X X DOUBLE PRECISION z X X PhiSec = DLOG( z ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PsiSec( w ) X X DOUBLE PRECISION w X X PsiSec = DEXP( w ) X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION PhiPriRecSec( z ) X X DOUBLE PRECISION z X X PhiPriRecSec = z X X RETURN X END X* ----+----------------------------------------------------------------+ X DOUBLE PRECISION FUNCTION RhoSec( z ) X X DOUBLE PRECISION z X X RhoSec = z X X RETURN X END X* ----+----------------------------------------------------------------+ X SUBROUTINE SincMNRD(H, M, N, X, Work) X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* This section defines arguments. X INTEGER M, N X DOUBLE PRECISION H, X, Work(*) X X* This section defines local constants. X DOUBLE PRECISION Pi X PARAMETER ( Pi = 3.141592653589792D0 ) X X DOUBLE PRECISION a0 X PARAMETER( a0 = 1.0D0 ) X X DOUBLE PRECISION a2 X PARAMETER ( a2 = 0.1666666716337204D0 ) X X DOUBLE PRECISION a4 X PARAMETER ( a4 = 0.0083333337679505D0 ) X X DOUBLE PRECISION a6 X PARAMETER ( a6 = 0.0001984127011383D0 ) X X DOUBLE PRECISION a8 X PARAMETER ( a8 = 0.000002755731922D0 ) X X* This section defines local variables. X INTEGER i, ii X DOUBLE PRECISION ww, ww2, sww X X* This is the start of the executable code. X ww = (Pi/H) * (X - DFLOAT(M) * H) X sww = DSIN(ww) X X ii = 0 X DO i = M, N X ii = ii + 1 X X IF ( DABS( ww ) .LE. 0.1 ) THEN X ww2 = ww * ww X Work(SincPointer + ii) = X . (((a8 * ww2 - a6) * ww2 + a4) * ww2 - a2) * ww2 + a0 X ELSE X Work(SincPointer + ii) = sww / ww X END IF X X ww = ww - Pi X sww = -sww X X END DO X X RETURN X END X X INTEGER FUNCTION Newton(H, N, NumOfEqua, InitC, Work, X . TotalAdd, MaxIter ) X X* n e w t o n . f X X* This procedure implements the sinc indefinite integral algorithm. X* It solves to an initial-value problem (IVP) of the form: X X* / X* | X* y = Y0 + | f( x, y) dx X* | X* / X X* by applying Newton's method to the equation, using sinc X* collocation to approximate the integrals. This code X* handles a vector of simultaneous IVP's. Since any system of X* simultaneous first order differential equations (ODE's) whose X* constraints are specified as initial conditions can be expressed as X* indefinite integral equations and since any higher order ODE may X* be represented as a system of first order ODE's, this software can X* handle a vary large class of problems. X X* The accuracy of a solution is governed by the number of "sinc points" X* used in approximating the solution. in the code below the number X* of sinc points = 2*N+1. The expected accuracy of the solution is: X X* Tol = EXP(-C*DSQRT(N)) X X* where Tol is a bound on absolute error between the approximated X* solution and the actual solution. Practical values for N are X* typically from 20 to 50, etc. Using N = 20, typically One should X* expect to 3 to 5 digits of accuracy (i.e., the absolute error should X* be on the order of 10^(3) to 10^(5) ). Using N = 50, One should X* expect to obtain 5 to 8 digits of accuracy. The program has a X* current limit: N < 201. X X* ARGUMENT LIST: X X* H - DOUBLE PRECISION - This is the step size. In some instances X* picking the optimal H can make a significant difference in the X* accuracy of the result, particularly for small N. However, as X* long H is picked as: X X* H = C / sqrt(N) for some C X X* the algorithm will have the appropriate asymptotic accuracy. X X* N - INTEGER - The size of the sinc approximation. The total X* number of terms in the approximation is actually 2N + 1. X X* NumOfEqua - INTEGER - Number of dependent variables in the system X* (i.e., the number of simultaneous IVP's to be solved). X X* InitC(NumOfEqua) - DOUBLE PRECISION array - The Initial values of X* the dependent variables. X X* Work - DOUBLE PRECISION array - Containing all `dynamically' X* allocated intermediate results. X X* TotalAdd - INTEGER - Indicates the amount of space used in Work. X* Probably this is a vestigial parameter ... but I'm making small X* steps. X X* MaxIter - INTEGER - The maximum number of iterations the procedure X* will use in trying to converge. X X* RESULT X X* IF convergence is achieved during the function returns 0. X* Otherwise the result is -1 and the values contained in the work X* array represent the approximation of the last iteration (typically X* of no value). X X* CREDITS: X X* This code was written by Michael O'Reilly, Brian Keyes, and X* Kenneth Parker, at the University of Utah, March 20, 1991. It was X* rewritten by Kenneth Parker, Michael O'Reilly, and Brian Keyes, at X* the University of Utah, July 30, 1992. All work was performed X* under the direction of Frank Stenger X X X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* This section defines the argument list. X INTEGER N, NumOfEqua, TotalAdd, MaxIter X DOUBLE PRECISION H, InitC(*), Work(*) X X* This section defines local constants. X INTEGER UpDateRate, FortranOrder X DOUBLE PRECISION Zero, Pi X X PARAMETER ( UpDateRate = 1 ) X PARAMETER ( FortranOrder = 0 ) X X PARAMETER ( Zero = 0.D0 ) X PARAMETER ( Pi = 3.1415926535897932D0 ) X X* This section defines local variables. X INTEGER i, J, k, EqNum, IdCount, L, Temp2 X INTEGER FunCount, Error, UpDateCount X INTEGER OldYPointer, FPointer, YPointer X X DOUBLE PRECISION Tol, ErrorNorm, Temp1 X X* Start of executable code. X L = 2*N+1 X X* Set the expected absolute Tolerance of the computed solution. X* This number may not be set arbitrarily; it is a function of the X* value of N. X print*,' Number of Sinc-points ', N X Tol = DSQRT(DBLE(N))*DEXP(-Pi*DSQRT(DBLE(N)/2.0D0)) X X WRITE(*,'(/A/)') ' Newton iteration method' X WRITE(*,'(A,G11.3/)') ' Tolerance = ',Tol X X* Use the starting values for the solution. X DO k = 1, NumOfEqua X Work(Y0Pointer + k) = InitC(k) X END DO X X* Start Newton's iteration. X UpDateCount = 0 X IdCount = 0 X FunCount = 0 X ErrorNorm = 2 * Tol X DO WHILE ((IdCount .LT. MaxIter) .AND. (ErrorNorm .GT. Tol)) X X IdCount = IdCount + 1 X PRINT *,' Iteration #: ', IdCount-1 X X* Call the function. X DO i = 1, L X CALL F(NumOfEqua, i, Work) X END DO X FunCount = FunCount + L X X X* Update the jacobian of F if needed. X* And, of course at the same time, update A and the LU decomposition X* of A. X IF (UpDateCount .GT. 0) THEN X UpDateCount = UpDateCount - 1 X ELSE X UpDateCount = UpDateRate - 1 X DO EqNum = 1,L X CALL J_ij(NumOfEqua, EqNum, Work ) X END DO X CALL UpdateA(NumOfEqua, L, H, Work) X CALL dgefa X . (Work(APointer+1), NumOfEqua*(2*N+1), NumOfEqua*(2*N+1), X . Work(APivotPointer+1), Error) X END IF X X* Switch the pointers between the old and new y-values. X Temp2 = YPointerPointer X YPointerPointer = OldYPointerPointer X OldYPointerPointer = Temp2 X X* Update the y-values X* Y_i = Y0_i + h \sum_j \sigma_{j-i} (f_j/\phi_j) X DO k=1,NumOfEqua X YPointer = NINT(Work(YPointerPointer + k)) X FPointer = NINT(Work(FPointerPointer + k)) X X DO i=1,L X Work(YPointer+i) = Zero X X DO J = 1,L X Work(YPointer + i) = Work(YPointer + i) + X . Work(OffSigma + J-i + L) * X . Work(OneOverPhiPrimePointer + J) * Work(FPointer + J) X END DO X X Work(YPointer + i) = Work(Y0Pointer + k) + X . H*Work(YPointer + i) X X END DO X X END DO X X* Set B_i = Y_i - OldY_i X DO k = 1, NumOfEqua X YPointer = NINT(Work(YPointerPointer + k)) X OldYPointer = NINT(Work(OldYPointerPointer + k)) X X DO i = 1,L X Work(BPointer + (k-1)*L+i) = X . Work(YPointer + i) - Work(OldYPointer + i) X END DO X END DO X X* Compute the step (this is the x which satisfies Big_A x = b); X* store the result in b. X CALL dgesl X . (Work(APointer + 1), NumOfEqua*(2*N+1), NumOfEqua*(2*N+1), X . Work(APivotPointer + 1), Work(BPointer + 1), FortranOrder) X X* Update y, i.e., Y = OldY + B. X DO k = 1, NumOfEqua X YPointer = NINT(Work(YPointerPointer + k)) X OldYPointer = NINT(Work(OldYPointerPointer + k)) X X DO i = 1,L X Work(YPointer + i) = X . Work(BPointer + (k-1)*L + i) + Work(OldYPointer + i) X END DO X END DO X X* Check for convergence using the infinity norm, and report the X* error norm. X ErrorNorm = Zero X DO k=1,NumOfEqua X YPointer = NINT(Work(YPointerPointer + k)) X OldYPointer = NINT(Work(OldYPointerPointer + k)) X X DO i = 1,L X Temp1 = DABS(Work(YPointer + i) - Work(OldYPointer + i)) X IF (Temp1 .GT. ErrorNorm) THEN X ErrorNorm = Temp1 X END IF X END DO X END DO X X PRINT 12, ErrorNorm X 12 format(' ErrorNorm= ',D12.3) X* End Newton iteration. X END DO X X* Check to see if procedure converged. X IF (ErrorNorm .LT. Tol) THEN X PRINT *,' The Newton method converged in ', IdCount-1, X +' iterations' X PRINT *,' using ', FunCount, ' function evaluations' X PRINT *,' ' X Newton = 0 X ELSE X PRINT *,' Newton iteration failed to converge for' X PRINT *,' the current value of N within ', IdCount X PRINT *,' iterations.' X Newton = -1 X END IF X X RETURN X END X X X X* ----+----------------------------------------------------------------+ X X* UpdateA X X* This subroutine computes A ... I - h I^{(-1)} D( 1 / \phi') Bij. X X SUBROUTINE UpdateA(NumOfEqua, L, H, Work) X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* This section defines the argument list X INTEGER NumOfEqua, L X DOUBLE PRECISION H, Work(*) X X* This section defines local variables. X INTEGER i, j, k, SubColumn, SubRow, Pos, BijPointer X X* This is the executable code. X Pos = 0 X DO i = 1, NumOfEqua X X DO SubColumn = 1, L X X k = NumOfEqua * (i - 1) X DO j = 1, NumOfEqua X X k = k + 1 X BijPointer = NINT(Work(BijPointerPointer + k)) X X DO SubRow = 1, L X Pos = Pos + 1 X X Work(APointer + Pos) = X . - Work(OffSigma + SubColumn - SubRow + L)*H* X . Work(OneOverPhiPrimePointer + SubColumn) * X . Work(BijPointer + SubColumn) X X IF ((i .EQ. j) .AND. (SubRow .EQ. SubColumn)) THEN X Work(APointer + Pos) = Work(APointer + Pos) + 1 X END IF X X END DO X END DO X X END DO X END DO X X RETURN X END X X* ----+----------------------------------------------------------------+ X X* ----+----------------------------------------------------------------+ X X* Init2Work X X* Initializes the work array with values needed for the interpolation. X X SUBROUTINE Init2Work(Nxs, X, Work, TotalAdd) X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* Functions called X DOUBLE PRECISION Phi X EXTERNAL Phi X DOUBLE PRECISION Rho X EXTERNAL Rho X X* This section defines the argument list. X INTEGER Nxs, TotalAdd X X DOUBLE PRECISION X(*), Work(*) X X* This section defines local variables X INTEGER i, ii, TADD X X* Start of the executable code. X X* Initialize the needed Heap. X PhiPointer = TotalAdd X TADD = TotalAdd + Nxs X RhoPointer = TADD X X ii = 0 X DO i=1,Nxs X ii = ii + 1 X Work(PhiPointer + ii) = Phi(X(ii)) X Work(RhoPointer + ii) = Rho(X(ii)) X END DO X X PRINT *,' ' X PRINT *,' Loading completed.' X X RETURN X END X X X* ----+----------------------------------------------------------------+ X X* Init1Work X X* Initializes the work array with anything to be used in Newton. X* By keeping intermediate results in a work array, the temptation to X* recompute these results is minimized. X X* NOTE: This routine only uses Phi and OneOverPhiPrime from the X* problem specification. X X SUBROUTINE Init1Work(N, NN, H, Work, TotalAdd) X IMPLICIT NONE X X INCLUDE 'workspace.h' X X* Functions called X DOUBLE PRECISION Psi X EXTERNAL Psi X DOUBLE PRECISION OneOverPhiPrime X EXTERNAL OneOverPhiPrime X X* This section defines the argument list. X INTEGER N, NN, TotalAdd X X DOUBLE PRECISION Work(*) X DOUBLE PRECISION H X X* This section defines local constants. X DOUBLE PRECISION Half X PARAMETER ( Half = 0.5D0 ) X X DOUBLE PRECISION Sigma(201) X X INTEGER i X X DATA (SIGMA(I+1),I=0,37) / XC X . 0.0000000000000000E+00, 0.5894898722360839 , X . 0.4514116667901405 , 0.5330932376182722 , X . 0.4749696698836553 , 0.5201071641913087 , X . 0.4832052174977473 , 0.5144159971233055 , X . 0.4873742250578202 , 0.5112301526369977 , X . 0.4898881711538789 , 0.5091957420082168 , X . 0.4915683516686011 , 0.5077846578125663 , X . 0.4927702093748033 , 0.5067486944720118 , X . 0.4936724151782212 , 0.5059559079171770 , X . 0.4943745528336628 , 0.5053297104401591 , X . 0.4949364995706957 , 0.5048226073041721 , X . 0.4953964150848574 , 0.5044035851980506 , X . 0.4957797661366434 , 0.5040515358439449 , X . 0.4961041974885161 , 0.5037515950319094 , X . 0.4963823201653523 , 0.5034929932779406 , X . 0.4966233866311803 , 0.5032677369477720 , X . 0.4968343388552230 , 0.5030697682021877 , X . 0.4970204870276673 , 0.5028944125557015 , X . 0.4971859623359650 , 0.5027380053824386 / X data (sigma(i+1),i=38,75) / X . 0.4973340269269068 , 0.5025976332159279 , X . 0.4974672909775609 , 0.5024709506908217 , X . 0.4975878678043933 , 0.5023560485249016 , X . 0.4976974867059375 , 0.5022513566773183 , X . 0.4977975763910226 , 0.5021555722141233 , X . 0.4978893275646686 , 0.5020676048274522 , X . 0.4979737405030820 , 0.5019865351657797 , X . 0.4980516616563340 , 0.5019115825934454 , X . 0.4981238121215731 , 0.5018420799807533 , X . 0.4981908100179338 , 0.5017774537988254 , X . 0.4982531882343348 , 0.5017172082611159 , X . 0.4983114086292758 , 0.5016609125833092 , X . 0.4983658734834106 , 0.5016081906689419 , X . 0.4984169348055853 , 0.5015587126985118 , X . 0.4984649019474762 , 0.5015121882244937 , X . 0.4985100478749211 , 0.5014683604668267 , X . 0.4985526143645061 , 0.5014270015722505 , X . 0.4985928163343139 , 0.5013879086527001 , X . 0.4986308454725853 , 0.5013509004573851 / X data (sigma(i+1),i=76,95) / X . 0.4986668732935894 , 0.5013158145633707 , X . 0.4987010537234771 , 0.5012825049928007 , X . 0.4987335252983494 , 0.5012508401830422 , X . 0.4987644130407324 , 0.5012207012502380 , X . 0.4987938300680404 , 0.5011919804979528 , X . 0.4988218789766499 , 0.5011645801314821 , X . 0.4988486530372723 , 0.5011384111454787 , X . 0.4988742372309752 , 0.5011133923582401 , X . 0.4988987091500929 , 0.5010894495705807 , X . 0.4989221397841480 , 0.5010665148309350 / X X DATA (SIGMA(I+1),I=96,133) / X . 0.4989445942075478 , 0.5010445257913609 , X . 0.4989661321830883 , 0.5010234251415910 , X . 0.4989868086930457 , 0.5010031601103150 , X . 0.4990066744077986 , 0.5009836820245585 , X . 0.4990257761003836 , 0.5009649459194093 , X . 0.4990441570141277 , 0.5009469101915090 , X . 0.4990618571894377 , 0.5009295362906861 , X . 0.4990789137549436 , 0.5009127884449231 , X . 0.4990953611874528 , 0.5008966334145260 , X . 0.4991112315445431 , 0.5008810402719440 , X . 0.4991265546730973 , 0.5008659802041686 , X . 0.4991413583966329 , 0.5008514263350565 , X . 0.4991556686839016 , 0.5008373535652727 , X . 0.4991695098009059 , 0.5008237384278461 , X . 0.4991829044482064 , 0.5008105589575910 , X . 0.4991958738851533 , 0.5007977945728658 , X . 0.4992084380424717 , 0.5007854259683299 , X . 0.4992206156244573 , 0.5007734350175222 , X . 0.4992324242018795 , 0.5007618046842314 / X data (sigma(i+1),i=134,171) / X . 0.4992438802965679 , 0.5007505189417415 , X . 0.4992549994585343 , 0.5007395626991506 , X . 0.4992657963363902 , 0.5007289217340515 , X . 0.4992762847417280 , 0.5007185826309395 , X . 0.4992864777080626 , 0.5007085327247896 , X . 0.4992963875448610 , 0.5006987600493020 , X . 0.4993060258871316 , 0.5006892532893728 , X . 0.4993154037409938 , 0.5006800017373920 , X . 0.4993245315256009 , 0.5006709952530141 , X . 0.4993334191117545 , 0.5006622242260842 , X . 0.4993420758575095 , 0.5006536795424343 , X . 0.4993505106410414 , 0.5006453525522925 , X . 0.4993587318910163 , 0.5006372350410781 , X . 0.4993667476146842 , 0.5006293192023715 , X . 0.4993745654238911 , 0.5006215976128758 , X . 0.4993821925591876 , 0.5006140632091974 , X . 0.4993896359121953 , 0.5006067092662964 , X . 0.4993969020463755 , 0.5005995293774664 , X . 0.4994039972163319 , 0.5005925174357194 / X data (sigma(i+1),i=172,191) / X . 0.4994109273857661 , 0.5005856676164615 , X . 0.4994176982441967 , 0.5005789743613571 , X . 0.4994243152225357 , 0.5005724323632864 , X . 0.4994307835076177 , 0.5005660365523119 , X . 0.4994371080557585 , 0.5005597820825741 , X . 0.4994432936054206 , 0.5005536643200484 , X . 0.4994493446890528 , 0.5005476788310940 , X . 0.4994552656441658 , 0.5005418213717391 , X . 0.4994610606237005 , 0.5005360878776455 , X . 0.4994667336057422 , 0.5005304744547047 / X X data (sigma(i+1),i=192,200) / X . 0.4994722884026277D0 , 0.5005249773702170D0 , X . 0.4994777286694879D0 , 0.5005195930446155D0 , X . 0.4994830579122687D0 , 0.5005143180436920D0 , X . 0.4994882794952658D0 , 0.5005091490712937D0 , X . 0.4994933966482060D0 / XC X* This section defines local variables. X INTEGER ii, k, kk, m, l, lV, Add X X DOUBLE PRECISION KH X Xc print*,' N: ', N Xc print*,' NN: ', NN Xc print*,' H: ', H X* Start of executable code. X l = 2*N + 1 X lV = 2*l - 1 X Add = 2*N + 1 X X* Set up the work common block. X TotalAdd = 0 X PsiPointer = TotalAdd X TotalAdd = TotalAdd + Add X X OneOverPhiPrimePointer = TotalAdd X TotalAdd = TotalAdd + Add X X OffExpKH = TotalAdd X TotalAdd = TotalAdd + Add X X SincPointer = TotalAdd X TotalAdd = TotalAdd + Add X X YPointerPointer = TotalAdd Xc print*,' YPointerPointer: ', YPointerPointer X TotalAdd = TotalAdd + NN X X OldYPointerPointer = TotalAdd X TotalAdd = TotalAdd + NN X X FPointerPointer = TotalAdd X TotalAdd = TotalAdd + NN X X BijPointerPointer = TotalAdd X TotalAdd = TotalAdd + NN * NN X X APointer = TotalAdd X TotalAdd = TotalAdd + NN*NN*l*l X X BPointer = TotalAdd X TotalAdd = TotalAdd + NN * Add X X APivotPointer = TotalAdd X TotalAdd = TotalAdd + NN * Add X X* Fill in the pointer lists that are stored within the work array. X DO k = 1, NN X Work(YPointerPointer + k) = DBLE(TotalAdd) X TotalAdd = TotalAdd + Add X END DO X X DO k = 1, NN X Work(OldYPointerPointer + k) = DBLE(TotalAdd) X TotalAdd = TotalAdd + Add X END DO X X DO k = 1, NN X Work(FPointerPointer + k) = DBLE(TotalAdd) X TotalAdd = TotalAdd + Add X END DO X X kk = 0 X DO m = 1, NN X DO k = 1, NN X kk = kk + 1 X Work(BijPointerPointer + kk) = DBLE(TotalAdd) X TotalAdd = TotalAdd + Add X END DO X END DO X X Y0Pointer = TotalAdd X TotalAdd = TotalAdd + NN X X OffSigma = TotalAdd X TotalAdd = TotalAdd + 2*l-1 X X* Initialize values in the work common block. X* Compute exp(k*H). X ii = 0 X DO i = -N, N X ii = ii + 1 X KH = DFLOAT(i) * H X Work( OffExpKH + ii ) = DEXP( KH ) X END DO X X* Compute i^{(-1)} and store in the work array. X Work(OffSigma + l) = Half X DO i = l+1, lV X Work(OffSigma + i) = Half - Sigma(i - l + 1) X END DO X X DO i = l-1, 1, -1 X Work(OffSigma + i) = Half + Sigma(l - i + 1) X END DO X X* Load the appropriate values into the work array. This section X* only uses Psi and PhiPriRec for the list of sinc functions. X PRINT *,' ' X PRINT *,' Loading values for the domain.' X X ii = 0 X DO i = -N, N X ii = ii + 1 X KH = DFLOAT(i) * H X Work(PsiPointer + ii) = Psi(KH) X Work(OneOverPhiPrimePointer + ii) = X . OneOverPhiPrime(Work(PsiPointer + ii)) X END DO X X PRINT *,' ' X PRINT *,' Loading completed.' Xc print*,' YPointerPointer: ', YPointerPointer X X RETURN X END X X* ----+----------------------------------------------------------------+ X X SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) X INTEGER LDA,N,IPVT(1),INFO X DOUBLE PRECISION A(LDA,1) XC XC DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. XC XC DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED XC DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. XC (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) . XC XC ON ENTRY XC XC A DOUBLE PRECISION(LDA, N) XC THE MATRIX TO BE FACTORED. XC XC LDA INTEGER XC THE LEADING DIMENSION OF THE ARRAY A . XC XC N INTEGER XC THE ORDER OF THE MATRIX A . XC XC ON RETURN XC XC A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS XC WHICH WERE USED TO OBTAIN IT. XC THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE XC L IS A PRODUCT OF PERMUTATION AND UNIT LOWER XC TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. XC XC IPVT INTEGER(N) XC AN INTEGER VECTOR OF PIVOT INDICES. XC XC INFO INTEGER XC = 0 NORMAL VALUE. XC = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR XC CONDITION FOR THIS SUBROUTINE, BUT IT DOES XC INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO XC IF CALLED. USE RCOND IN DGECO FOR A RELIABLE XC INDICATION OF SINGULARITY. XC XC LINPACK. THIS VERSION DATED 08/14/78 . XC CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. XC XC SUBROUTINES AND FUNCTIONS XC XC BLAS DAXPY,DSCAL,IDAMAX XC XC INTERNAL VARIABLES XC X DOUBLE PRECISION T X INTEGER IDAMAX,J,K,KP1,L,NM1 XC XC XC GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING XC X INFO = 0 X NM1 = N - 1 X IF (NM1 .LT. 1) GO TO 70 X DO 60 K = 1, NM1 X KP1 = K + 1 XC XC FIND L = PIVOT INDEX XC X L = IDAMAX(N-K+1,A(K,K),1) + K - 1 X IPVT(K) = L XC XC ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED XC X IF (A(L,K) .EQ. 0.0D0) GO TO 40 XC XC INTERCHANGE IF NECESSARY XC X IF (L .EQ. K) GO TO 10 X T = A(L,K) X A(L,K) = A(K,K) X A(K,K) = T X 10 CONTINUE XC XC COMPUTE MULTIPLIERS XC X T = -1.0D0/A(K,K) X CALL DSCAL(N-K,T,A(K+1,K),1) XC XC ROW ELIMINATION WITH COLUMN INDEXING XC X DO 30 J = KP1, N X T = A(L,J) X IF (L .EQ. K) GO TO 20 X A(L,J) = A(K,J) X A(K,J) = T X 20 CONTINUE X CALL DAXPY(N-K,T,A(K+1,K),1,A(K+1,J),1) X 30 CONTINUE X GO TO 50 X 40 CONTINUE X INFO = K X 50 CONTINUE X 60 CONTINUE X 70 CONTINUE X IPVT(N) = N X IF (A(N,N) .EQ. 0.0D0) INFO = N X RETURN X END X X SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) X INTEGER LDA,N,IPVT(1),JOB X DOUBLE PRECISION A(LDA,1),B(1) XC XC DGESL SOLVES THE DOUBLE PRECISION SYSTEM XC A * X = B OR TRANS(A) * X = B XC USING THE FACTORS COMPUTED BY DGECO OR DGEFA. XC XC ON ENTRY XC XC A DOUBLE PRECISION(LDA, N) XC THE OUTPUT FROM DGECO OR DGEFA. XC XC LDA INTEGER XC THE LEADING DIMENSION OF THE ARRAY A . XC XC N INTEGER XC THE ORDER OF THE MATRIX A . XC XC IPVT INTEGER(N) XC THE PIVOT VECTOR FROM DGECO OR DGEFA. XC XC B DOUBLE PRECISION(N) XC THE RIGHT HAND SIDE VECTOR. XC XC JOB INTEGER XC = 0 TO SOLVE A*X = B , XC = NONZERO TO SOLVE TRANS(A)*X = B WHERE XC TRANS(A) IS THE TRANSPOSE. XC XC ON RETURN XC XC B THE SOLUTION VECTOR X . XC XC ERROR CONDITION XC XC A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A XC ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY XC BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER XC SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE XC CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 XC OR DGEFA HAS SET INFO .EQ. 0 . XC XC TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX XC WITH P COLUMNS XC CALL DGECO(A,LDA,N,IPVT,RCOND,Z) XC IF (RCOND IS TOO SMALL) GO TO ... XC DO 10 J = 1, P XC CALL DGESL(A,LDA,N,IPVT,C(1,J),0) XC 10 CONTINUE XC XC LINPACK. THIS VERSION DATED 08/14/78 . XC CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. XC XC SUBROUTINES AND FUNCTIONS XC XC BLAS DAXPY,DDOT XC XC INTERNAL VARIABLES XC X DOUBLE PRECISION DDOT,T X INTEGER K,KB,L,NM1 XC X NM1 = N - 1 X IF (JOB .NE. 0) GO TO 50 XC XC JOB = 0 , SOLVE A * X = B XC FIRST SOLVE L*Y = B XC X IF (NM1 .LT. 1) GO TO 30 X DO 20 K = 1, NM1 X L = IPVT(K) X T = B(L) X IF (L .EQ. K) GO TO 10 X B(L) = B(K) X B(K) = T X 10 CONTINUE X CALL DAXPY(N-K,T,A(K+1,K),1,B(K+1),1) X 20 CONTINUE X 30 CONTINUE XC XC NOW SOLVE U*X = Y XC X DO 40 KB = 1, N X K = N + 1 - KB X B(K) = B(K)/A(K,K) X T = -B(K) X CALL DAXPY(K-1,T,A(1,K),1,B(1),1) X 40 CONTINUE X GO TO 100 X 50 CONTINUE XC XC JOB = NONZERO, SOLVE TRANS(A) * X = B XC FIRST SOLVE TRANS(U)*Y = B XC X DO 60 K = 1, N X T = DDOT(K-1,A(1,K),1,B(1),1) X B(K) = (B(K) - T)/A(K,K) X 60 CONTINUE XC XC NOW SOLVE TRANS(L)*X = Y XC X IF (NM1 .LT. 1) GO TO 90 X DO 80 KB = 1, NM1 X K = N - KB X B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) X L = IPVT(K) X IF (L .EQ. K) GO TO 70 X T = B(L) X B(L) = B(K) X B(K) = T X 70 CONTINUE X 80 CONTINUE X 90 CONTINUE X 100 CONTINUE X RETURN X END X DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) XC XC FORMS THE DOT PRODUCT OF TWO VECTORS. XC USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. XC JACK DONGARRA, LINPACK, 3/11/78. XC X DOUBLE PRECISION DX(1),DY(1),DTEMP X INTEGER I,INCX,INCY,IX,IY,M,MP1,N XC X DDOT = 0.0D0 X DTEMP = 0.0D0 X IF(N.LE.0)RETURN X IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 XC XC CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS XC NOT EQUAL TO 1 XC X IX = 1 X IY = 1 X IF(INCX.LT.0)IX = (-N+1)*INCX + 1 X IF(INCY.LT.0)IY = (-N+1)*INCY + 1 X DO 10 I = 1,N X DTEMP = DTEMP + DX(IX)*DY(IY) X IX = IX + INCX X IY = IY + INCY X 10 CONTINUE X DDOT = DTEMP X RETURN XC XC CODE FOR BOTH INCREMENTS EQUAL TO 1 XC XC XC CLEAN-UP LOOP XC X 20 M = MOD(N,5) X IF( M .EQ. 0 ) GO TO 40 X DO 30 I = 1,M X DTEMP = DTEMP + DX(I)*DY(I) X 30 CONTINUE X IF( N .LT. 5 ) GO TO 60 X 40 MP1 = M + 1 X DO 50 I = MP1,N,5 X DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + X * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) X 50 CONTINUE X 60 DDOT = DTEMP X RETURN X END X SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) XC XC CONSTANT TIMES A VECTOR PLUS A VECTOR. XC USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. XC JACK DONGARRA, LINPACK, 3/11/78. XC X DOUBLE PRECISION DX(1),DY(1),DA X INTEGER I,IX,IY,INCX,INCY,M,MP1,N XC X IF(N.LE.0)RETURN X IF (DA .EQ. 0.0D0) RETURN X IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 XC XC CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS XC NOT EQUAL TO 1 XC X IX = 1 X IY = 1 X IF(INCX.LT.0)IX = (-N+1)*INCX + 1 X IF(INCY.LT.0)IY = (-N+1)*INCY + 1 X DO 10 I = 1,N X DY(IY) = DY(IY) + DA*DX(IX) X IX = IX + INCX X IY = IY + INCY X 10 CONTINUE X RETURN XC XC CODE FOR BOTH INCREMENTS EQUAL TO 1 XC XC XC CLEAN-UP LOOP XC X 20 M = MOD(N,4) X IF( M .EQ. 0 ) GO TO 40 X DO 30 I = 1,M X DY(I) = DY(I) + DA*DX(I) X 30 CONTINUE X IF( N .LT. 4 ) RETURN X 40 MP1 = M + 1 X DO 50 I = MP1,N,4 X DY(I) = DY(I) + DA*DX(I) X DY(I + 1) = DY(I + 1) + DA*DX(I + 1) X DY(I + 2) = DY(I + 2) + DA*DX(I + 2) X DY(I + 3) = DY(I + 3) + DA*DX(I + 3) X 50 CONTINUE X RETURN X END X SUBROUTINE DSCAL(N,DA,DX,INCX) XC XC SCALES A VECTOR BY A CONSTANT. XC USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. XC JACK DONGARRA, LINPACK, 3/11/78. XC X DOUBLE PRECISION DA,DX(1) X INTEGER I,INCX,M,MP1,N,NINCX XC X IF(N.LE.0)RETURN X IF(INCX.EQ.1)GO TO 20 XC XC CODE FOR INCREMENT NOT EQUAL TO 1 XC X NINCX = N*INCX X DO 10 I = 1,NINCX,INCX X DX(I) = DA*DX(I) X 10 CONTINUE X RETURN XC XC CODE FOR INCREMENT EQUAL TO 1 XC XC XC CLEAN-UP LOOP XC X 20 M = MOD(N,5) X IF( M .EQ. 0 ) GO TO 40 X DO 30 I = 1,M X DX(I) = DA*DX(I) X 30 CONTINUE X IF( N .LT. 5 ) RETURN X 40 MP1 = M + 1 X DO 50 I = MP1,N,5 X DX(I) = DA*DX(I) X DX(I + 1) = DA*DX(I + 1) X DX(I + 2) = DA*DX(I + 2) X DX(I + 3) = DA*DX(I + 3) X DX(I + 4) = DA*DX(I + 4) X 50 CONTINUE X RETURN X END X INTEGER FUNCTION IDAMAX(N,DX,INCX) XC XC FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. XC JACK DONGARRA, LINPACK, 3/11/78. XC X DOUBLE PRECISION DX(1),DMAX X INTEGER I,INCX,IX,N XC X IDAMAX = 0 X IF( N .LT. 1 ) RETURN X IDAMAX = 1 X IF(N.EQ.1)RETURN X IF(INCX.EQ.1)GO TO 20 XC XC CODE FOR INCREMENT NOT EQUAL TO 1 XC X IX = 1 X DMAX = DABS(DX(1)) X IX = IX + INCX X DO 10 I = 2,N X IF(DABS(DX(IX)).LE.DMAX) GO TO 5 X IDAMAX = I X DMAX = DABS(DX(IX)) X 5 IX = IX + INCX X 10 CONTINUE X RETURN XC XC CODE FOR INCREMENT EQUAL TO 1 XC X 20 DMAX = DABS(DX(1)) X DO 30 I = 2,N X IF(DABS(DX(I)).LE.DMAX) GO TO 30 X IDAMAX = I X DMAX = DABS(DX(I)) X 30 CONTINUE X RETURN X END X X SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) X INTEGER LDA,N,IPVT(1) X DOUBLE PRECISION A(LDA,1),Z(1) X DOUBLE PRECISION RCOND XC XC DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION XC AND ESTIMATES THE CONDITION OF THE MATRIX. XC XC IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. XC TO SOLVE A*X = B , FOLLOW DGECO BY DGESL. XC TO COMPUTE INVERSE(A)*C , FOLLOW DGECO BY DGESL. XC TO COMPUTE DETERMINANT(A) , FOLLOW DGECO BY DGEDI. XC TO COMPUTE INVERSE(A) , FOLLOW DGECO BY DGEDI. XC XC ON ENTRY XC XC A DOUBLE PRECISION(LDA, N) XC THE MATRIX TO BE FACTORED. XC XC LDA INTEGER XC THE LEADING DIMENSION OF THE ARRAY A . XC XC N INTEGER XC THE ORDER OF THE MATRIX A . XC XC ON RETURN XC XC A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS XC WHICH WERE USED TO OBTAIN IT. XC THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE XC L IS A PRODUCT OF PERMUTATION AND UNIT LOWER XC TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. XC XC IPVT INTEGER(N) XC AN INTEGER VECTOR OF PIVOT INDICES. XC XC RCOND DOUBLE PRECISION XC AN ESTIMATE OF THE RECIPROCAL CONDITION OF A . XC FOR THE SYSTEM A*X = B , RELATIVE PERTURBATIONS XC IN A AND B OF SIZE EPSILON MAY CAUSE XC RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND . XC IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION XC 1.0 + RCOND .EQ. 1.0 XC IS TRUE, THEN A MAY BE SINGULAR TO WORKING XC PRECISION. IN PARTICULAR, RCOND IS ZERO IF XC EXACT SINGULARITY IS DETECTED OR THE ESTIMATE XC UNDERFLOWS. XC XC Z DOUBLE PRECISION(N) XC A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. XC IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS XC AN APPROXIMATE NULL VECTOR IN THE SENSE THAT XC NORM(A*Z) = RCOND*NORM(A)*NORM(Z) . XC XC LINPACK. THIS VERSION DATED 08/14/78 . XC CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB. XC XC SUBROUTINES AND FUNCTIONS XC XC LINPACK DGEFA XC BLAS DAXPY,DDOT,DSCAL,DASUM XC FORTRAN DABS,DMAX1,DSIGN XC XC INTERNAL VARIABLES XC X DOUBLE PRECISION DDOT,EK,T,WK,WKM X DOUBLE PRECISION ANORM,S,DASUM,SM,YNORM X INTEGER INFO,J,K,KB,KP1,L XC XC XC COMPUTE 1-NORM OF A XC X ANORM = 0.0D0 X DO 10 J = 1, N X ANORM = DMAX1(ANORM,DASUM(N,A(1,J),1)) X 10 CONTINUE XC XC FACTOR XC X CALL DGEFA(A,LDA,N,IPVT,INFO) XC XC RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) . XC ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E . XC TRANS(A) IS THE TRANSPOSE OF A . THE COMPONENTS OF E ARE XC CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE XC TRANS(U)*W = E . THE VECTORS ARE FREQUENTLY RESCALED TO AVOID XC OVERFLOW. XC XC SOLVE TRANS(U)*W = E XC X EK = 1.0D0 X DO 20 J = 1, N X Z(J) = 0.0D0 X 20 CONTINUE X DO 100 K = 1, N X IF (Z(K) .NE. 0.0D0) EK = DSIGN(EK,-Z(K)) X IF (DABS(EK-Z(K)) .LE. DABS(A(K,K))) GO TO 30 X S = DABS(A(K,K))/DABS(EK-Z(K)) X CALL DSCAL(N,S,Z,1) X EK = S*EK X 30 CONTINUE X WK = EK - Z(K) X WKM = -EK - Z(K) X S = DABS(WK) X SM = DABS(WKM) X IF (A(K,K) .EQ. 0.0D0) GO TO 40 X WK = WK/A(K,K) X WKM = WKM/A(K,K) X GO TO 50 X 40 CONTINUE X WK = 1.0D0 X WKM = 1.0D0 X 50 CONTINUE X KP1 = K + 1 X IF (KP1 .GT. N) GO TO 90 X DO 60 J = KP1, N X SM = SM + DABS(Z(J)+WKM*A(K,J)) X Z(J) = Z(J) + WK*A(K,J) X S = S + DABS(Z(J)) X 60 CONTINUE X IF (S .GE. SM) GO TO 80 X T = WKM - WK X WK = WKM X DO 70 J = KP1, N X Z(J) = Z(J) + T*A(K,J) X 70 CONTINUE X 80 CONTINUE X 90 CONTINUE X Z(K) = WK X 100 CONTINUE X S = 1.0D0/DASUM(N,Z,1) X CALL DSCAL(N,S,Z,1) XC XC SOLVE TRANS(L)*Y = W XC X DO 120 KB = 1, N X K = N + 1 - KB X IF (K .LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) X IF (DABS(Z(K)) .LE. 1.0D0) GO TO 110 X S = 1.0D0/DABS(Z(K)) X CALL DSCAL(N,S,Z,1) X 110 CONTINUE X L = IPVT(K) X T = Z(L) X Z(L) = Z(K) X Z(K) = T X 120 CONTINUE X S = 1.0D0/DASUM(N,Z,1) X CALL DSCAL(N,S,Z,1) XC X YNORM = 1.0D0 XC XC SOLVE L*V = Y XC X DO 140 K = 1, N X L = IPVT(K) X T = Z(L) X Z(L) = Z(K) X Z(K) = T X IF (K .LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+1),1) X IF (DABS(Z(K)) .LE. 1.0D0) GO TO 130 X S = 1.0D0/DABS(Z(K)) X CALL DSCAL(N,S,Z,1) X YNORM = S*YNORM X 130 CONTINUE X 140 CONTINUE X S = 1.0D0/DASUM(N,Z,1) X CALL DSCAL(N,S,Z,1) X YNORM = S*YNORM XC XC SOLVE U*Z = V XC X DO 160 KB = 1, N X K = N + 1 - KB X IF (DABS(Z(K)) .LE. DABS(A(K,K))) GO TO 150 X S = DABS(A(K,K))/DABS(Z(K)) X CALL DSCAL(N,S,Z,1) X YNORM = S*YNORM X 150 CONTINUE X IF (A(K,K) .NE. 0.0D0) Z(K) = Z(K)/A(K,K) X IF (A(K,K) .EQ. 0.0D0) Z(K) = 1.0D0 X T = -Z(K) X CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) X 160 CONTINUE XC MAKE ZNORM = 1.0 X S = 1.0D0/DASUM(N,Z,1) X CALL DSCAL(N,S,Z,1) X YNORM = S*YNORM XC X IF (ANORM .NE. 0.0D0) RCOND = YNORM/ANORM X IF (ANORM .EQ. 0.0D0) RCOND = 0.0D0 X RETURN X END X DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) XC XC TAKES THE SUM OF THE ABSOLUTE VALUES. XC JACK DONGARRA, LINPACK, 3/11/78. XC X DOUBLE PRECISION DX(1),DTEMP X INTEGER I,INCX,M,MP1,N,NINCX XC X DASUM = 0.0D0 X DTEMP = 0.0D0 X IF(N.LE.0)RETURN X IF(INCX.EQ.1)GO TO 20 XC XC CODE FOR INCREMENT NOT EQUAL TO 1 XC X NINCX = N*INCX X DO 10 I = 1,NINCX,INCX X DTEMP = DTEMP + DABS(DX(I)) X 10 CONTINUE X DASUM = DTEMP X RETURN XC XC CODE FOR INCREMENT EQUAL TO 1 XC XC XC CLEAN-UP LOOP XC X 20 M = MOD(N,6) X IF( M .EQ. 0 ) GO TO 40 X DO 30 I = 1,M X DTEMP = DTEMP + DABS(DX(I)) X 30 CONTINUE X IF( N .LT. 6 ) GO TO 60 X 40 MP1 = M + 1 X DO 50 I = MP1,N,6 X DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) X * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) X 50 CONTINUE X 60 DASUM = DTEMP X RETURN X END X SUBROUTINE CALJY0(ARG,RESULT,JINT) XC--------------------------------------------------------------------- XC XC This packet computes zero-order Bessel functions of the first and XC second kind (J0 and Y0), for real arguments X, where 0 < X <= XMAX XC for Y0, and |X| <= XMAX for J0. It contains two function-type XC subprograms, BESJ0 and BESY0, and one subroutine-type XC subprogram, CALJY0. The calling statements for the primary XC entries are: XC XC Y = BESJ0(X) XC and XC Y = BESY0(X), XC XC where the entry points correspond to the functions J0(X) and Y0(X), XC respectively. The routine CALJY0 is intended for internal packet XC use only, all computations within the packet being concentrated in XC this one routine. The function subprograms invoke CALJY0 with XC the statement XC CALL CALJY0(ARG,RESULT,JINT), XC where the parameter usage is as follows: XC XC Function Parameters for CALJY0 XC call ARG RESULT JINT XC XC BESJ0(ARG) |ARG| .LE. XMAX J0(ARG) 0 XC BESY0(ARG) 0 .LT. ARG .LE. XMAX Y0(ARG) 1 XC XC The main computation uses unpublished minimax rational XC approximations for X .LE. 8.0, and an approximation from the XC book Computer Approximations by Hart, et. al., Wiley and Sons, XC New York, 1968, for arguments larger than 8.0 Part of this XC transportable packet is patterned after the machine-dependent XC FUNPACK program BESJ0(X), but cannot match that version for XC efficiency or accuracy. This version uses rational functions XC that are theoretically accurate to at least 18 significant decimal XC digits for X <= 8, and at least 18 decimal places for X > 8. The XC accuracy achieved depends on the arithmetic system, the compiler, XC the intrinsic functions, and proper selection of the machine- XC dependent constants. XC XC******************************************************************* XC XC Explanation of machine-dependent constants XC XC XINF = largest positive machine number XC XMAX = largest acceptable argument. The functions AINT, SIN XC and COS must perform properly for ABS(X) .LE. XMAX. XC We recommend that XMAX be a small integer multiple of XC sqrt(1/eps), where eps is the smallest positive number XC such that 1+eps > 1. XC XSMALL = positive argument such that 1.0-(X/2)**2 = 1.0 XC to machine precision for all ABS(X) .LE. XSMALL. XC We recommend that XSMALL < sqrt(eps)/beta, where beta XC is the floating-point radix (usually 2 or 16). XC XC Approximate values for some important machines are XC XC eps XMAX XSMALL XINF XC XC CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322 XC CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465 XC IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38 XC IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308 XC IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75 XC UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307 XC VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38 XC XC******************************************************************* XC******************************************************************* XC XC Error Returns XC XC The program returns the value zero for X .GT. XMAX, and returns XC -XINF when BESLY0 is called with a negative or zero argument. XC XC XC Intrinsic functions required are: XC XC ABS, AINT, COS, LOG, SIN, SQRT XC XC XC Latest modification: June 2, 1989 XC XC Author: W. J. Cody XC Mathematics and Computer Science Division XC Argonne National Laboratory XC Argonne, IL 60439 XC XC-------------------------------------------------------------------- X INTEGER I,JINT XCS REAL X DOUBLE PRECISION X 1 ARG,AX,CONS,DOWN,EIGHT,FIVE5,FOUR,ONE,ONEOV8,PI2,PJ0, X 2 PJ1,PLG,PROD,PY0,PY1,PY2,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1, X 3 QY2,Q0,Q1,RESJ,RESULT,R0,R1,SIXTY4,THREE,TWOPI,TWOPI1, X 4 TWOPI2,TWO56,UP,W,WSQ,XDEN,XINF,XMAX,XNUM,XSMALL,XJ0, X 5 XJ1,XJ01,XJ02,XJ11,XJ12,XY,XY0,XY01,XY02,XY1,XY11,XY12, X 6 XY2,XY21,XY22,Z,ZERO,ZSQ X DIMENSION PJ0(7),PJ1(8),PLG(4),PY0(6),PY1(7),PY2(8),P0(6),P1(6), X 1 QJ0(5),QJ1(7),QLG(4),QY0(5),QY1(6),QY2(7),Q0(5),Q1(5) XC------------------------------------------------------------------- XC Mathematical constants XC CONS = ln(.5) + Euler's gamma XC------------------------------------------------------------------- XCS DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0E0,1.0E0,3.0E0,4.0E0,8.0E0/, XCS 1 FIVE5,SIXTY4,ONEOV8,P17/5.5E0,64.0E0,0.125E0,1.716E-1/, XCS 2 TWO56,CONS/256.0E0,-1.1593151565841244881E-1/, XCS 3 PI2,TWOPI/6.3661977236758134308E-1,6.2831853071795864769E0/, XCS 4 TWOPI1,TWOPI2/6.28125E0,1.9353071795864769253E-3/ X DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0D0,1.0D0,3.0D0,4.0D0,8.0D0/, X 1 FIVE5,SIXTY4,ONEOV8,P17/5.5D0,64.0D0,0.125D0,1.716D-1/, X 2 TWO56,CONS/256.0D0,-1.1593151565841244881D-1/, X 3 PI2,TWOPI/6.3661977236758134308D-1,6.2831853071795864769D0/, X 4 TWOPI1,TWOPI2/6.28125D0,1.9353071795864769253D-3/ XC------------------------------------------------------------------- XC Machine-dependent constants XC------------------------------------------------------------------- XCS DATA XMAX/8.19E+03/,XSMALL/1.22E-09/,XINF/1.7E+38/ X DATA XMAX/1.07D+09/,XSMALL/9.31D-10/,XINF/1.7D+38/ XC------------------------------------------------------------------- XC Zeroes of Bessel functions XC------------------------------------------------------------------- XCS DATA XJ0/2.4048255576957727686E+0/,XJ1/5.5200781102863106496E+0/, XCS 1 XY0/8.9357696627916752158E-1/,XY1/3.9576784193148578684E+0/, XCS 2 XY2/7.0860510603017726976E+0/, XCS 3 XJ01/ 616.0E+0/, XJ02/-1.4244423042272313784E-03/, XCS 4 XJ11/1413.0E+0/, XJ12/ 5.4686028631064959660E-04/, XCS 5 XY01/ 228.0E+0/, XY02/ 2.9519662791675215849E-03/, XCS 6 XY11/1013.0E+0/, XY12/ 6.4716931485786837568E-04/, XCS 7 XY21/1814.0E+0/, XY22/ 1.1356030177269762362E-04/ X DATA XJ0/2.4048255576957727686D+0/,XJ1/5.5200781102863106496D+0/, X 1 XY0/8.9357696627916752158D-1/,XY1/3.9576784193148578684D+0/, X 2 XY2/7.0860510603017726976D+0/, X 3 XJ01/ 616.0D+0/, XJ02/-1.4244423042272313784D-03/, X 4 XJ11/1413.0D+0/, XJ12/ 5.4686028631064959660D-04/, X 5 XY01/ 228.0D+0/, XY02/ 2.9519662791675215849D-03/, X 6 XY11/1013.0D+0/, XY12/ 6.4716931485786837568D-04/, X 7 XY21/1814.0D+0/, XY22/ 1.1356030177269762362D-04/ XC------------------------------------------------------------------- XC Coefficients for rational approximation to ln(x/a) XC-------------------------------------------------------------------- XCS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02, XCS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/ XCS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02, XCS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/ X DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02, X 1 -5.4989956895857911039D+02,3.5687548468071500413D+02/ X DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02, X 1 -3.3442903192607538956D+02,1.7843774234035750207D+02/ XC------------------------------------------------------------------- XC Coefficients for rational approximation of XC J0(X) / (X**2 - XJ0**2), XSMALL < |X| <= 4.0 XC-------------------------------------------------------------------- XCS DATA PJ0/6.6302997904833794242E+06,-6.2140700423540120665E+08, XCS 1 2.7282507878605942706E+10,-4.1298668500990866786E+11, XCS 2 -1.2117036164593528341E-01, 1.0344222815443188943E+02, XCS 3 -3.6629814655107086448E+04/ XCS DATA QJ0/4.5612696224219938200E+05, 1.3985097372263433271E+08, XCS 1 2.6328198300859648632E+10, 2.3883787996332290397E+12, XCS 2 9.3614022392337710626E+02/ X DATA PJ0/6.6302997904833794242D+06,-6.2140700423540120665D+08, X 1 2.7282507878605942706D+10,-4.1298668500990866786D+11, X 2 -1.2117036164593528341D-01, 1.0344222815443188943D+02, X 3 -3.6629814655107086448D+04/ X DATA QJ0/4.5612696224219938200D+05, 1.3985097372263433271D+08, X 1 2.6328198300859648632D+10, 2.3883787996332290397D+12, X 2 9.3614022392337710626D+02/ XC------------------------------------------------------------------- XC Coefficients for rational approximation of XC J0(X) / (X**2 - XJ1**2), 4.0 < |X| <= 8.0 XC------------------------------------------------------------------- XCS DATA PJ1/4.4176707025325087628E+03, 1.1725046279757103576E+04, XCS 1 1.0341910641583726701E+04,-7.2879702464464618998E+03, XCS 2 -1.2254078161378989535E+04,-1.8319397969392084011E+03, XCS 3 4.8591703355916499363E+01, 7.4321196680624245801E+02/ XCS DATA QJ1/3.3307310774649071172E+02,-2.9458766545509337327E+03, XCS 1 1.8680990008359188352E+04,-8.4055062591169562211E+04, XCS 2 2.4599102262586308984E+05,-3.5783478026152301072E+05, XCS 3 -2.5258076240801555057E+01/ X DATA PJ1/4.4176707025325087628D+03, 1.1725046279757103576D+04, X 1 1.0341910641583726701D+04,-7.2879702464464618998D+03, X 2 -1.2254078161378989535D+04,-1.8319397969392084011D+03, X 3 4.8591703355916499363D+01, 7.4321196680624245801D+02/ X DATA QJ1/3.3307310774649071172D+02,-2.9458766545509337327D+03, X 1 1.8680990008359188352D+04,-8.4055062591169562211D+04, X 2 2.4599102262586308984D+05,-3.5783478026152301072D+05, X 3 -2.5258076240801555057D+01/ XC------------------------------------------------------------------- XC Coefficients for rational approximation of XC (Y0(X) - 2 LN(X/XY0) J0(X)) / (X**2 - XY0**2), XC XSMALL < |X| <= 3.0 XC-------------------------------------------------------------------- XCS DATA PY0/1.0102532948020907590E+04,-2.1287548474401797963E+06, XCS 1 2.0422274357376619816E+08,-8.3716255451260504098E+09, XCS 2 1.0723538782003176831E+11,-1.8402381979244993524E+01/ XCS DATA QY0/6.6475986689240190091E+02, 2.3889393209447253406E+05, XCS 1 5.5662956624278251596E+07, 8.1617187777290363573E+09, XCS 2 5.8873865738997033405E+11/ X DATA PY0/1.0102532948020907590D+04,-2.1287548474401797963D+06, X 1 2.0422274357376619816D+08,-8.3716255451260504098D+09, X 2 1.0723538782003176831D+11,-1.8402381979244993524D+01/ X DATA QY0/6.6475986689240190091D+02, 2.3889393209447253406D+05, X 1 5.5662956624278251596D+07, 8.1617187777290363573D+09, X 2 5.8873865738997033405D+11/ XC------------------------------------------------------------------- XC Coefficients for rational approximation of XC (Y0(X) - 2 LN(X/XY1) J0(X)) / (X**2 - XY1**2), XC 3.0 < |X| <= 5.5 XC-------------------------------------------------------------------- XCS DATA PY1/-1.4566865832663635920E+04, 4.6905288611678631510E+06, XCS 1 -6.9590439394619619534E+08, 4.3600098638603061642E+10, XCS 2 -5.5107435206722644429E+11,-2.2213976967566192242E+13, XCS 3 1.7427031242901594547E+01/ XCS DATA QY1/ 8.3030857612070288823E+02, 4.0669982352539552018E+05, XCS 1 1.3960202770986831075E+08, 3.4015103849971240096E+10, XCS 2 5.4266824419412347550E+12, 4.3386146580707264428E+14/ X DATA PY1/-1.4566865832663635920D+04, 4.6905288611678631510D+06, X 1 -6.9590439394619619534D+08, 4.3600098638603061642D+10, X 2 -5.5107435206722644429D+11,-2.2213976967566192242D+13, X 3 1.7427031242901594547D+01/ X DATA QY1/ 8.3030857612070288823D+02, 4.0669982352539552018D+05, X 1 1.3960202770986831075D+08, 3.4015103849971240096D+10, X 2 5.4266824419412347550D+12, 4.3386146580707264428D+14/ XC------------------------------------------------------------------- XC Coefficients for rational approximation of XC (Y0(X) - 2 LN(X/XY2) J0(X)) / (X**2 - XY2**2), XC 5.5 < |X| <= 8.0 XC-------------------------------------------------------------------- XCS DATA PY2/ 2.1363534169313901632E+04,-1.0085539923498211426E+07, XCS 1 2.1958827170518100757E+09,-1.9363051266772083678E+11, XCS 2 -1.2829912364088687306E+11, 6.7016641869173237784E+14, XCS 3 -8.0728726905150210443E+15,-1.7439661319197499338E+01/ XCS DATA QY2/ 8.7903362168128450017E+02, 5.3924739209768057030E+05, XCS 1 2.4727219475672302327E+08, 8.6926121104209825246E+10, XCS 2 2.2598377924042897629E+13, 3.9272425569640309819E+15, XCS 3 3.4563724628846457519E+17/ X DATA PY2/ 2.1363534169313901632D+04,-1.0085539923498211426D+07, X 1 2.1958827170518100757D+09,-1.9363051266772083678D+11, X 2 -1.2829912364088687306D+11, 6.7016641869173237784D+14, X 3 -8.0728726905150210443D+15,-1.7439661319197499338D+01/ X DATA QY2/ 8.7903362168128450017D+02, 5.3924739209768057030D+05, X 1 2.4727219475672302327D+08, 8.6926121104209825246D+10, X 2 2.2598377924042897629D+13, 3.9272425569640309819D+15, X 3 3.4563724628846457519D+17/ XC------------------------------------------------------------------- XC Coefficients for Hart,s approximation, |X| > 8.0 XC------------------------------------------------------------------- XCS DATA P0/3.4806486443249270347E+03, 2.1170523380864944322E+04, XCS 1 4.1345386639580765797E+04, 2.2779090197304684302E+04, XCS 2 8.8961548424210455236E-01, 1.5376201909008354296E+02/ XCS DATA Q0/3.5028735138235608207E+03, 2.1215350561880115730E+04, XCS 1 4.1370412495510416640E+04, 2.2779090197304684318E+04, XCS 2 1.5711159858080893649E+02/ XCS DATA P1/-2.2300261666214198472E+01,-1.1183429920482737611E+02, XCS 1 -1.8591953644342993800E+02,-8.9226600200800094098E+01, XCS 2 -8.8033303048680751817E-03,-1.2441026745835638459E+00/ XCS DATA Q1/1.4887231232283756582E+03, 7.2642780169211018836E+03, XCS 1 1.1951131543434613647E+04, 5.7105024128512061905E+03, XCS 2 9.0593769594993125859E+01/ X DATA P0/3.4806486443249270347D+03, 2.1170523380864944322D+04, X 1 4.1345386639580765797D+04, 2.2779090197304684302D+04, X 2 8.8961548424210455236D-01, 1.5376201909008354296D+02/ X DATA Q0/3.5028735138235608207D+03, 2.1215350561880115730D+04, X 1 4.1370412495510416640D+04, 2.2779090197304684318D+04, X 2 1.5711159858080893649D+02/ X DATA P1/-2.2300261666214198472D+01,-1.1183429920482737611D+02, X 1 -1.8591953644342993800D+02,-8.9226600200800094098D+01, X 2 -8.8033303048680751817D-03,-1.2441026745835638459D+00/ X DATA Q1/1.4887231232283756582D+03, 7.2642780169211018836D+03, X 1 1.1951131543434613647D+04, 5.7105024128512061905D+03, X 2 9.0593769594993125859D+01/ XC------------------------------------------------------------------- XC Check for error conditions XC------------------------------------------------------------------- X AX = ABS(ARG) X IF ((JINT .EQ. 1) .AND. (ARG .LE. ZERO)) THEN X RESULT = -XINF X GO TO 2000 X ELSE IF (AX .GT. XMAX) THEN X RESULT = ZERO X GO TO 2000 X END IF X IF (AX .GT. EIGHT) GO TO 800 X IF (AX .LE. XSMALL) THEN X IF (JINT .EQ. 0) THEN X RESULT = ONE X ELSE X RESULT = PI2 * (LOG(AX) + CONS) X END IF X GO TO 2000 X END IF XC------------------------------------------------------------------- XC Calculate J0 for appropriate interval, preserving XC accuracy near the zero of J0 XC------------------------------------------------------------------- X ZSQ = AX * AX X IF (AX .LE. FOUR) THEN X XNUM = (PJ0(5) * ZSQ + PJ0(6)) * ZSQ + PJ0(7) X XDEN = ZSQ + QJ0(5) X DO 50 I = 1, 4 X XNUM = XNUM * ZSQ + PJ0(I) X XDEN = XDEN * ZSQ + QJ0(I) X 50 CONTINUE X PROD = ((AX - XJ01/TWO56) - XJ02) * (AX + XJ0) X ELSE X WSQ = ONE - ZSQ / SIXTY4 X XNUM = PJ1(7) * WSQ + PJ1(8) X XDEN = WSQ + QJ1(7) X DO 220 I = 1, 6 X XNUM = XNUM * WSQ + PJ1(I) X XDEN = XDEN * WSQ + QJ1(I) X 220 CONTINUE X PROD = (AX + XJ1) * ((AX - XJ11/TWO56) - XJ12) X END IF X RESULT = PROD * XNUM / XDEN X IF (JINT .EQ. 0) GO TO 2000 XC------------------------------------------------------------------- XC Calculate Y0. First find RESJ = pi/2 ln(x/xn) J0(x), XC where xn is a zero of Y0 XC------------------------------------------------------------------- X IF (AX .LE. THREE) THEN X UP = (AX-XY01/TWO56)-XY02 X XY = XY0 X ELSE IF (AX .LE. FIVE5) THEN X UP = (AX-XY11/TWO56)-XY12 X XY = XY1 X ELSE X UP = (AX-XY21/TWO56)-XY22 X XY = XY2 X END IF X DOWN = AX + XY X IF (ABS(UP) .LT. P17*DOWN) THEN X W = UP/DOWN X WSQ = W*W X XNUM = PLG(1) X XDEN = WSQ + QLG(1) X DO 320 I = 2, 4 X XNUM = XNUM*WSQ + PLG(I) X XDEN = XDEN*WSQ + QLG(I) X 320 CONTINUE X RESJ = PI2 * RESULT * W * XNUM/XDEN X ELSE X RESJ = PI2 * RESULT * LOG(AX/XY) X END IF XC------------------------------------------------------------------- XC Now calculate Y0 for appropriate interval, preserving XC accuracy near the zero of Y0 XC------------------------------------------------------------------- X IF (AX .LE. THREE) THEN X XNUM = PY0(6) * ZSQ + PY0(1) X XDEN = ZSQ + QY0(1) X DO 340 I = 2, 5 X XNUM = XNUM * ZSQ + PY0(I) X XDEN = XDEN * ZSQ + QY0(I) X 340 CONTINUE X ELSE IF (AX .LE. FIVE5) THEN X XNUM = PY1(7) * ZSQ + PY1(1) X XDEN = ZSQ + QY1(1) X DO 360 I = 2, 6 X XNUM = XNUM * ZSQ + PY1(I) X XDEN = XDEN * ZSQ + QY1(I) X 360 CONTINUE X ELSE X XNUM = PY2(8) * ZSQ + PY2(1) X XDEN = ZSQ + QY2(1) X DO 380 I = 2, 7 X XNUM = XNUM * ZSQ + PY2(I) X XDEN = XDEN * ZSQ + QY2(I) X 380 CONTINUE X END IF X RESULT = RESJ + UP * DOWN * XNUM / XDEN X GO TO 2000 XC------------------------------------------------------------------- XC Calculate J0 or Y0 for |ARG| > 8.0 XC------------------------------------------------------------------- X 800 Z = EIGHT / AX X W = AX / TWOPI X W = AINT(W) + ONEOV8 X W = (AX - W * TWOPI1) - W * TWOPI2 X ZSQ = Z * Z X XNUM = P0(5) * ZSQ + P0(6) X XDEN = ZSQ + Q0(5) X UP = P1(5) * ZSQ + P1(6) X DOWN = ZSQ + Q1(5) X DO 850 I = 1, 4 X XNUM = XNUM * ZSQ + P0(I) X XDEN = XDEN * ZSQ + Q0(I) X UP = UP * ZSQ + P1(I) X DOWN = DOWN * ZSQ + Q1(I) X 850 CONTINUE X R0 = XNUM / XDEN X R1 = UP / DOWN X IF (JINT .EQ. 0) THEN X RESULT = SQRT(PI2/AX) * (R0*COS(W) - Z*R1*SIN(W)) X ELSE X RESULT = SQRT(PI2/AX) * (R0*SIN(W) + Z*R1*COS(W)) X END IF X 2000 RETURN XC---------- Last line of CALJY0 ---------- X END X X DOUBLE PRECISION FUNCTION BESJ0(X) XCS REAL FUNCTION BESJ0(X) XC-------------------------------------------------------------------- XC XC This subprogram computes approximate values for Bessel functions XC of the first kind of order zero for arguments |X| <= XMAX XC (see comments heading CALJY0). XC XC-------------------------------------------------------------------- X INTEGER JINT XCS REAL X, RESULT X DOUBLE PRECISION X, RESULT XC-------------------------------------------------------------------- X JINT=0 X CALL CALJY0(X,RESULT,JINT) X BESJ0 = RESULT X RETURN XC---------- Last line of BESJ0 ---------- X END X DOUBLE PRECISION FUNCTION BESY0(X) XCS REAL FUNCTION BESY0(X) XC-------------------------------------------------------------------- XC XC This subprogram computes approximate values for Bessel functions XC of the second kind of order zero for arguments 0 < X <= XMAX XC (see comments heading CALJY0). XC XC-------------------------------------------------------------------- X INTEGER JINT XCS REAL X, RESULT X DOUBLE PRECISION X, RESULT XC-------------------------------------------------------------------- X JINT=1 X CALL CALJY0(X,RESULT,JINT) X BESY0 = RESULT X RETURN XC---------- Last line of BESY0 ---------- X END X X* ====================================================================== X* NIST Guide to Available Math Software. X* Fullsource for module DGAMI from package CMLIB. X* Retrieved from CAMSUN on Tue May 23 18:39:46 1995. X* ====================================================================== X DOUBLE PRECISION FUNCTION DGAMI(A,X) XC***BEGIN PROLOGUE DGAMI XC***DATE WRITTEN 770701 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. C7E XC***KEYWORDS DOUBLE PRECISION,GAMMA,GAMMA FUNCTION, XC INCOMPLETE GAMMA FUNCTION,SPECIAL FUNCTION XC***AUTHOR FULLERTON, W., (LANL) XC***PURPOSE Calculates the d.p. incomplete Gamma function. XC***DESCRIPTION XC XC Evaluate the incomplete gamma function defined by XC XC DGAMI = integral from T = 0 to X of EXP(-T) * T**(A-1.0) . XC XC DGAMI is evaluated for positive values of A and non-negative values XC of X. A slight deterioration of 2 or 3 digits accuracy will occur XC when DGAMI is very large or very small, because logarithmic variables XC are used. The function and both arguments are double precision. XC***REFERENCES (NONE) XC***ROUTINES CALLED DGAMIT,DLNGAM,XERROR XC***END PROLOGUE DGAMI X DOUBLE PRECISION A, X, FACTOR, DLNGAM, DGAMIT XC***FIRST EXECUTABLE STATEMENT DGAMI X IF (A.LE.0.D0) CALL XERROR ( 'DGAMI A MUST BE GT ZERO', 25, 1,2) X IF (X.LT.0.D0) CALL XERROR ( 'DGAMI X MUST BE GE ZERO', 25, 2,2) XC X DGAMI = 0.D0 X IF (X.EQ.0.0D0) RETURN XC XC THE ONLY ERROR POSSIBLE IN THE EXPRESSION BELOW IS A FATAL OVERFLOW. X FACTOR = DEXP (DLNGAM(A) + A*DLOG(X)) XC X DGAMI = FACTOR * DGAMIT (A, X) XC X RETURN X END X DOUBLE PRECISION FUNCTION DGAMIT(A,X) XC***BEGIN PROLOGUE DGAMIT XC***DATE WRITTEN 770701 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. C7E XC***KEYWORDS COMPLEMENTARY,COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, XC DOUBLE PRECISION,GAMMA FUNCTION,SPECIAL FUNCTION,TRICOMI XC***AUTHOR FULLERTON, W., (LANL) XC***PURPOSE Calculates Tricomi's form of the incomplete Gamma function. XC***DESCRIPTION XC XC Evaluate Tricomi's incomplete gamma function defined by XC XC DGAMIT = X**(-A)/GAMMA(A) * integral T = 0 to X of EXP(-T) * T**(A-1.) XC XC for A .GT. 0.0 and by analytic XC continuation for A .LE. 0.0. Gamma(X) is the complete XC gamma function of X. DGAMIT is evaluated for arbitrary real values of XC A and for non-negative values of X (even though DGAMIT is defined for XC X .LT. 0.0), except that for X = 0 and A .LE. 0.0, DGAMIT is infinite, XC a fatal error. The function and both arguments are double precision. XC XC A slight deterioration of 2 or 3 digits accuracy will occur when XC DGAMIT is very large or very small in absolute value, because log- XC arithmic variables are used. Also, if the parameter A is very close XC to a negative integer (but not a negative integer), there is a loss XC of accuracy, which is reported if the result is less than half XC machine precision. XC XC Ref. -- W. Gautschi, An Evaluation Procedure for Incomplete Gamma XC Functions, ACM Trans. Math. Software, Vol. 5, No. 4, December 1979. XC***REFERENCES (NONE) XC***ROUTINES CALLED D1MACH,D9GMIT,D9LGIC,D9LGIT,DGAMR,DINT,DLGAMS, XC DLNGAM,XERCLR,XERROR XC***END PROLOGUE DGAMIT X DOUBLE PRECISION A, X, AEPS, AINTA, ALGAP1, ALNEPS, ALNG, ALX, X 1 BOT, H, SGA, SGNGAM, SQEPS, T, D1MACH, DGAMR, D9GMIT, D9LGIT, X 2 DLNGAM, D9LGIC, DINT X DATA ALNEPS, SQEPS, BOT / 3*0.D0 / XC***FIRST EXECUTABLE STATEMENT DGAMIT X IF (ALNEPS.NE.0.D0) GO TO 10 X ALNEPS = -DLOG (D1MACH(3)) X SQEPS = DSQRT (D1MACH(4)) X BOT = DLOG (D1MACH(1)) XC X 10 IF (X.LT.0.D0) CALL XERROR ( 'DGAMIT X IS NEGATIVE', 21, 2, 2) XC X IF (X.NE.0.D0) ALX = DLOG (X) X SGA = 1.0D0 X IF (A.NE.0.D0) SGA = DSIGN (1.0D0, A) X AINTA = DINT (A + 0.5D0*SGA) X AEPS = A - AINTA XC X IF (X.GT.0.D0) GO TO 20 X DGAMIT = 0.0D0 X IF (AINTA.GT.0.D0 .OR. AEPS.NE.0.D0) DGAMIT = DGAMR(A+1.0D0) X RETURN XC X 20 IF (X.GT.1.D0) GO TO 30 X IF (A.GE.(-0.5D0) .OR. AEPS.NE.0.D0) CALL DLGAMS (A+1.0D0, ALGAP1, X 1 SGNGAM) X DGAMIT = D9GMIT (A, X, ALGAP1, SGNGAM, ALX) X RETURN XC X 30 IF (A.LT.X) GO TO 40 X T = D9LGIT (A, X, DLNGAM(A+1.0D0)) X IF (T.LT.BOT) CALL XERCLR X DGAMIT = DEXP (T) X RETURN XC X 40 ALNG = D9LGIC (A, X, ALX) XC XC EVALUATE DGAMIT IN TERMS OF DLOG (DGAMIC (A, X)) XC X H = 1.0D0 X IF (AEPS.EQ.0.D0 .AND. AINTA.LE.0.D0) GO TO 50 XC X CALL DLGAMS (A+1.0D0, ALGAP1, SGNGAM) X T = DLOG (DABS(A)) + ALNG - ALGAP1 X IF (T.GT.ALNEPS) GO TO 60 XC X IF (T.GT.(-ALNEPS)) H = 1.0D0 - SGA * SGNGAM * DEXP(T) X IF (DABS(H).GT.SQEPS) GO TO 50 XC X CALL XERCLR X CALL XERROR ( 'DGAMIT RESULT LT HALF PRECISION', 32, 1, 1) XC X 50 T = -A*ALX + DLOG(DABS(H)) X IF (T.LT.BOT) CALL XERCLR X DGAMIT = DSIGN (DEXP(T), H) X RETURN XC X 60 T = T - A*ALX X IF (T.LT.BOT) CALL XERCLR X DGAMIT = -SGA * SGNGAM * DEXP(T) X RETURN XC X END X DOUBLE PRECISION FUNCTION DGAMR(X) XC***BEGIN PROLOGUE DGAMR XC***DATE WRITTEN 770701 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. C7A XC***KEYWORDS DOUBLE PRECISION,GAMMA,GAMMA FUNCTION, XC RECIPROCAL GAMMA FUNCTION,SPECIAL FUNCTION XC***AUTHOR FULLERTON, W., (LANL) XC***PURPOSE Calculates d.p. reciprocal Gamma function. XC***DESCRIPTION XC XC DGAMR(X) calculates the double precision reciprocal of the XC complete gamma function for double precision argument X. XC***REFERENCES (NONE) XC***ROUTINES CALLED DGAMMA,DINT,DLGAMS,XERCLR,XGETF,XSETF XC***END PROLOGUE DGAMR X EXTERNAL DGAMMA X INTEGER IROLD X DOUBLE PRECISION X, ALNGX, SGNGX, DGAMMA, DINT XC***FIRST EXECUTABLE STATEMENT DGAMR X DGAMR = 0.0D0 X IF (X.LE.0.0D0 .AND. DINT(X).EQ.X) RETURN XC X CALL XGETF (IROLD) X CALL XSETF (1) X IF (DABS(X).GT.10.0D0) GO TO 10 X DGAMR = 1.0D0/DGAMMA(X) X CALL XERCLR X CALL XSETF (IROLD) X RETURN XC X 10 CALL DLGAMS (X, ALNGX, SGNGX) X CALL XERCLR X CALL XSETF (IROLD) X DGAMR = SGNGX * DEXP(-ALNGX) X RETURN XC X END X SUBROUTINE DLGAMS(X,DLGAM,SGNGAM) XC***BEGIN PROLOGUE DLGAMS XC***DATE WRITTEN 770701 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. C7A XC***KEYWORDS ABSOLUTE VALUE,DOUBLE PRECISION,GAMMA FUNCTION,LOGARITHM, XC SPECIAL FUNCTION XC***AUTHOR FULLERTON, W., (LANL) XC***PURPOSE Calculates the log of the absolute value of the Gamma XC function XC***DESCRIPTION XC XC DLGAMS(X,DLGAM,SGNGAM) calculates the double precision natural XC logarithm of the absolute value of the gamma function for XC double precision argument X and stores the result in double XC precision argument DLGAM. XC***REFERENCES (NONE) XC***ROUTINES CALLED DINT,DLNGAM XC***END PROLOGUE DLGAMS X DOUBLE PRECISION X, DLGAM, SGNGAM, DINT, DLNGAM, INT XC***FIRST EXECUTABLE STATEMENT DLGAMS X DLGAM = DLNGAM(X) X SGNGAM = 1.0D0 X IF (X.GT.0.D0) RETURN XC X INT = DMOD (-DINT(X), 2.0D0) + 0.1D0 X IF (INT.EQ.0) SGNGAM = -1.0D0 XC X RETURN X END X DOUBLE PRECISION FUNCTION DLNGAM(X) XC***BEGIN PROLOGUE DLNGAM XC***DATE WRITTEN 770601 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. C7A XC***KEYWORDS ABSOLUTE VALUE,DOUBLE PRECISION,GAMMA FUNCTION,LOGARITHM, XC SPECIAL FUNCTION XC***AUTHOR FULLERTON, W., (LANL) XC***PURPOSE Computes the d.p. logarithm of the absolute value of the XC Gamma function XC***DESCRIPTION XC XC DLNGAM(X) calculates the double precision logarithm of the XC absolute value of the gamma function for double precision XC argument X. XC***REFERENCES (NONE) XC***ROUTINES CALLED D1MACH,D9LGMC,DGAMMA,DINT,XERROR XC***END PROLOGUE DLNGAM X EXTERNAL DGAMMA X DOUBLE PRECISION X, DXREL, PI, SINPIY, SQPI2L, SQ2PIL, XMAX, X 1 Y, DINT, DGAMMA, D9LGMC, D1MACH,TEMP X DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 / X DATA SQPI2L / +.2257913526 4472743236 3097614947 441 D+0 / X DATA PI / 3.1415926535 8979323846 2643383279 50 D0 / X DATA XMAX, DXREL / 2*0.D0 / XC***FIRST EXECUTABLE STATEMENT DLNGAM X IF (XMAX.NE.0.D0) GO TO 10 X TEMP = 1.0D0/DLOG(D1MACH(2)) X XMAX = TEMP * D1MACH(2) X DXREL = DSQRT (D1MACH(4)) XC X 10 Y = DABS (X) X IF (Y.GT.10.D0) GO TO 20 XC XC DLOG (DABS (DGAMMA(X)) ) FOR DABS(X) .LE. 10.0 XC X DLNGAM = DLOG (DABS (DGAMMA(X)) ) X RETURN XC XC DLOG ( DABS (DGAMMA(X)) ) FOR DABS(X) .GT. 10.0 XC X 20 IF (Y.GT.XMAX) CALL XERROR ( 'DLNGAM DABS(X) SO BIG DLNGAM OVERFL X 1OWS', 39, 2, 2) XC X IF (X.GT.0.D0) DLNGAM = SQ2PIL + (X-0.5D0)*DLOG(X) - X + D9LGMC(Y) X IF (X.GT.0.D0) RETURN XC X SINPIY = DABS (DSIN(PI*Y)) X IF (SINPIY.EQ.0.D0) CALL XERROR ( 'DLNGAM X IS A NEGATIVE INTEGER X 1', 31, 3, 2) XC X IF (DABS ((X-DINT(X-0.5D0))/X).LT.DXREL) CALL XERROR ( 'DLNGAM X 1ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER', X 2,68,1, 1) XC X DLNGAM = SQPI2L + (X-0.5D0)*DLOG(Y) - X - DLOG(SINPIY) - D9LGMC(Y) X RETURN XC X END X SUBROUTINE XERCLR XC***BEGIN PROLOGUE XERCLR XC***DATE WRITTEN 790801 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3C XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Resets current error number to zero. XC***DESCRIPTION XC Abstract XC This routine simply resets the current error number to zero. XC This may be necessary to do in order to determine that XC a certain error has occurred again since the last time XC NUMXER was referenced. XC XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC Latest revision --- 7 June 1978 XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED J4SAVE XC***END PROLOGUE XERCLR XC***FIRST EXECUTABLE STATEMENT XERCLR X INTEGER J4SAVE X INTEGER JUNK X JUNK = J4SAVE(1,0,.TRUE.) X RETURN X END X SUBROUTINE XERROR(MESSG,NMESSG,NERR,LEVEL) XC***BEGIN PROLOGUE XERROR XC***DATE WRITTEN 790801 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3C XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Processes an error (diagnostic) message. XC***DESCRIPTION XC Abstract XC XERROR processes a diagnostic message, in a manner XC determined by the value of LEVEL and the current value XC of the library error control flag, KONTRL. XC (See subroutine XSETF for details.) XC XC Description of Parameters XC --Input-- XC MESSG - the Hollerith message to be processed, containing XC no more than 72 characters. XC NMESSG- the actual number of characters in MESSG. XC NERR - the error number associated with this message. XC NERR must not be zero. XC LEVEL - error category. XC =2 means this is an unconditionally fatal error. XC =1 means this is a recoverable error. (I.e., it is XC non-fatal if XSETF has been appropriately called.) XC =0 means this is a warning message only. XC =-1 means this is a warning message which is to be XC printed at most once, regardless of how many XC times this call is executed. XC XC Examples XC CALL XERROR('SMOOTH -- NUM WAS ZERO.',23,1,2) XC CALL XERROR('INTEG -- LESS THAN FULL ACCURACY ACHIEVED.', XC 43,2,1) XC CALL XERROR('ROOTER -- ACTUAL ZERO OF F FOUND BEFORE INTERVAL XC 1FULLY COLLAPSED.',65,3,0) XC CALL XERROR('EXP -- UNDERFLOWS BEING SET TO ZERO.',39,1,-1) XC XC Latest revision --- 19 MAR 1980 XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED XERRWV XC***END PROLOGUE XERROR X CHARACTER*(*) MESSG X INTEGER LEVEL, NERR, NMESSG XC***FIRST EXECUTABLE STATEMENT XERROR X CALL XERRWV(MESSG,NMESSG,NERR,LEVEL,0,0,0,0,0.,0.) X RETURN X END X SUBROUTINE XERRWV(MESSG,NMESSG,NERR,LEVEL,NI,I1,I2,NR,R1,R2) XC***BEGIN PROLOGUE XERRWV XC***DATE WRITTEN 800319 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3C XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Processes error message allowing 2 integer and two real XC values to be included in the message. XC***DESCRIPTION XC Abstract XC XERRWV processes a diagnostic message, in a manner XC determined by the value of LEVEL and the current value XC of the library error control flag, KONTRL. XC (See subroutine XSETF for details.) XC In addition, up to two integer values and two real XC values may be printed along with the message. XC XC Description of Parameters XC --Input-- XC MESSG - the Hollerith message to be processed. XC NMESSG- the actual number of characters in MESSG. XC NERR - the error number associated with this message. XC NERR must not be zero. XC LEVEL - error category. XC =2 means this is an unconditionally fatal error. XC =1 means this is a recoverable error. (I.e., it is XC non-fatal if XSETF has been appropriately called.) XC =0 means this is a warning message only. XC =-1 means this is a warning message which is to be XC printed at most once, regardless of how many XC times this call is executed. XC NI - number of integer values to be printed. (0 to 2) XC I1 - first integer value. XC I2 - second integer value. XC NR - number of real values to be printed. (0 to 2) XC R1 - first real value. XC R2 - second real value. XC XC Examples XC CALL XERRWV('SMOOTH -- NUM (=I1) WAS ZERO.',29,1,2, XC 1 1,NUM,0,0,0.,0.) XC CALL XERRWV('QUADXY -- REQUESTED ERROR (R1) LESS THAN MINIMUM ( XC 1R2).,54,77,1,0,0,0,2,ERRREQ,ERRMIN) XC XC Latest revision --- 19 MAR 1980 XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED FDUMP,I1MACH,J4SAVE,XERABT,XERCTL,XERPRT,XERSAV, XC XGETUA XC***END PROLOGUE XERRWV X CHARACTER*(*) MESSG X CHARACTER*20 LFIRST X CHARACTER*37 FORM X INTEGER NMESSG, NERR, NI, I1,I2,NR, LUN X INTEGER J4SAVE, I1MACH X DIMENSION LUN(5) X DOUBLE PRECISION R1, R2 X INTEGER LKNTRL, MAXMES, KDUMMY, JUNK, KOUNT, LMESSG, LERR X INTEGER LLEVEL, MKNTRL,NUNIT, ISIZEI,ISIZEF, KUNIT,IUNIT X INTEGER I, IFATAL, LEVEL XC GET FLAGS XC***FIRST EXECUTABLE STATEMENT XERRWV X LKNTRL = J4SAVE(2,0,.FALSE.) X MAXMES = J4SAVE(4,0,.FALSE.) XC CHECK FOR VALID INPUT X IF ((NMESSG.GT.0).AND.(NERR.NE.0).AND. X 1 (LEVEL.GE.(-1)).AND.(LEVEL.LE.2)) GO TO 10 X IF (LKNTRL.GT.0) CALL XERPRT('FATAL ERROR IN...',17) X CALL XERPRT('XERROR -- INVALID INPUT',23) X IF (LKNTRL.GT.0) CALL FDUMP X IF (LKNTRL.GT.0) CALL XERPRT('JOB ABORT DUE TO FATAL ERROR.', X 1 29) X IF (LKNTRL.GT.0) CALL XERSAV(' ',0,0,0,KDUMMY) X CALL XERABT('XERROR -- INVALID INPUT',23) X RETURN X 10 CONTINUE XC RECORD MESSAGE X JUNK = J4SAVE(1,NERR,.TRUE.) X CALL XERSAV(MESSG,NMESSG,NERR,LEVEL,KOUNT) XC LET USER OVERRIDE X LFIRST = MESSG X LMESSG = NMESSG X LERR = NERR X LLEVEL = LEVEL X CALL XERCTL(LFIRST,LMESSG,LERR,LLEVEL,LKNTRL) XC RESET TO ORIGINAL VALUES X LMESSG = NMESSG X LERR = NERR X LLEVEL = LEVEL X LKNTRL = MAX0(-2,MIN0(2,LKNTRL)) X MKNTRL = IABS(LKNTRL) XC DECIDE WHETHER TO PRINT MESSAGE X IF ((LLEVEL.LT.2).AND.(LKNTRL.EQ.0)) GO TO 100 X IF (((LLEVEL.EQ.(-1)).AND.(KOUNT.GT.MIN0(1,MAXMES))) X 1.OR.((LLEVEL.EQ.0) .AND.(KOUNT.GT.MAXMES)) X 2.OR.((LLEVEL.EQ.1) .AND.(KOUNT.GT.MAXMES).AND.(MKNTRL.EQ.1)) X 3.OR.((LLEVEL.EQ.2) .AND.(KOUNT.GT.MAX0(1,MAXMES)))) GO TO 100 X IF (LKNTRL.LE.0) GO TO 20 X CALL XERPRT(' ',1) XC INTRODUCTION X IF (LLEVEL.EQ.(-1)) CALL XERPRT X 1('WARNING MESSAGE...THIS MESSAGE WILL ONLY BE PRINTED ONCE.',57) X IF (LLEVEL.EQ.0) CALL XERPRT('WARNING IN...',13) X IF (LLEVEL.EQ.1) CALL XERPRT X 1 ('RECOVERABLE ERROR IN...',23) X IF (LLEVEL.EQ.2) CALL XERPRT('FATAL ERROR IN...',17) X 20 CONTINUE XC MESSAGE X CALL XERPRT(MESSG,LMESSG) X CALL XGETUA(LUN,NUNIT) X ISIZEI = LOG10(FLOAT(I1MACH(9))) + 1.0 X ISIZEF = LOG10(FLOAT(I1MACH(10))**I1MACH(11)) + 1.0 X DO 50 KUNIT=1,NUNIT X IUNIT = LUN(KUNIT) X IF (IUNIT.EQ.0) IUNIT = I1MACH(4) X DO 22 I=1,MIN(NI,2) X WRITE (FORM,21) I,ISIZEI X 21 FORMAT ('(11X,21HIN ABOVE MESSAGE, I',I1,'=,I',I2,') ') X IF (I.EQ.1) WRITE (IUNIT,FORM) I1 X IF (I.EQ.2) WRITE (IUNIT,FORM) I2 X 22 CONTINUE X DO 24 I=1,MIN(NR,2) X WRITE (FORM,23) I,ISIZEF+10,ISIZEF X 23 FORMAT ('(11X,21HIN ABOVE MESSAGE, R',I1,'=,E', X 1 I2,'.',I2,')') X IF (I.EQ.1) WRITE (IUNIT,FORM) R1 X IF (I.EQ.2) WRITE (IUNIT,FORM) R2 X 24 CONTINUE X IF (LKNTRL.LE.0) GO TO 40 XC ERROR NUMBER X WRITE (IUNIT,30) LERR X 30 FORMAT (15H ERROR NUMBER =,I10) X 40 CONTINUE X 50 CONTINUE XC TRACE-BACK X IF (LKNTRL.GT.0) CALL FDUMP X 100 CONTINUE X IFATAL = 0 X IF ((LLEVEL.EQ.2).OR.((LLEVEL.EQ.1).AND.(MKNTRL.EQ.2))) X 1IFATAL = 1 XC QUIT HERE IF MESSAGE IS NOT FATAL X IF (IFATAL.LE.0) RETURN X IF ((LKNTRL.LE.0).OR.(KOUNT.GT.MAX0(1,MAXMES))) GO TO 120 XC PRINT REASON FOR ABORT X IF (LLEVEL.EQ.1) CALL XERPRT X 1 ('JOB ABORT DUE TO UNRECOVERED ERROR.',35) X IF (LLEVEL.EQ.2) CALL XERPRT X 1 ('JOB ABORT DUE TO FATAL ERROR.',29) XC PRINT ERROR SUMMARY X CALL XERSAV(' ',-1,0,0,KDUMMY) X 120 CONTINUE XC ABORT X IF ((LLEVEL.EQ.2).AND.(KOUNT.GT.MAX0(1,MAXMES))) LMESSG = 0 X CALL XERABT(MESSG,LMESSG) X RETURN X END X SUBROUTINE XERSAV(MESSG,NMESSG,NERR,LEVEL,ICOUNT) XC***BEGIN PROLOGUE XERSAV XC***DATE WRITTEN 800319 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. Z XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Records that an error occurred. XC***DESCRIPTION XC Abstract XC Record that this error occurred. XC XC Description of Parameters XC --Input-- XC MESSG, NMESSG, NERR, LEVEL are as in XERROR, XC except that when NMESSG=0 the tables will be XC dumped and cleared, and when NMESSG is less than zero the XC tables will be dumped and not cleared. XC --Output-- XC ICOUNT will be the number of times this message has XC been seen, or zero if the table has overflowed and XC does not contain this message specifically. XC When NMESSG=0, ICOUNT will not be altered. XC XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC Latest revision --- 19 Mar 1980 XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED I1MACH,XGETUA XC***END PROLOGUE XERSAV X INTEGER I1MACH X INTEGER I, II,ICOUNT, IUNIT, KUNIT, KOUNT, KOUNTX, LEVEL X INTEGER LEVTAB, NERR, NERTAB, NUNIT, NMESSG X INTEGER LUN(5) X CHARACTER*(*) MESSG X CHARACTER*20 MESTAB(10),MES X DIMENSION NERTAB(10),LEVTAB(10),KOUNT(10) X SAVE MESTAB,NERTAB,LEVTAB,KOUNT,KOUNTX XC NEXT TWO DATA STATEMENTS ARE NECESSARY TO PROVIDE A BLANK XC ERROR TABLE INITIALLY X DATA KOUNT(1),KOUNT(2),KOUNT(3),KOUNT(4),KOUNT(5), X 1 KOUNT(6),KOUNT(7),KOUNT(8),KOUNT(9),KOUNT(10) X 2 /0,0,0,0,0,0,0,0,0,0/ X DATA KOUNTX/0/ XC***FIRST EXECUTABLE STATEMENT XERSAV X IF (NMESSG.GT.0) GO TO 80 XC DUMP THE TABLE X IF (KOUNT(1).EQ.0) RETURN XC PRINT TO EACH UNIT X CALL XGETUA(LUN,NUNIT) X DO 60 KUNIT=1,NUNIT X IUNIT = LUN(KUNIT) X IF (IUNIT.EQ.0) IUNIT = I1MACH(4) XC PRINT TABLE HEADER X WRITE (IUNIT,10) X 10 FORMAT (32H0 ERROR MESSAGE SUMMARY/ X 1 51H MESSAGE START NERR LEVEL COUNT) XC PRINT BODY OF TABLE X DO 20 I=1,10 X IF (KOUNT(I).EQ.0) GO TO 30 X WRITE (IUNIT,15) MESTAB(I),NERTAB(I),LEVTAB(I),KOUNT(I) X 15 FORMAT (1X,A20,3I10) X 20 CONTINUE X 30 CONTINUE XC PRINT NUMBER OF OTHER ERRORS X IF (KOUNTX.NE.0) WRITE (IUNIT,40) KOUNTX X 40 FORMAT (41H0OTHER ERRORS NOT INDIVIDUALLY TABULATED=,I10) X WRITE (IUNIT,50) X 50 FORMAT (1X) X 60 CONTINUE X IF (NMESSG.LT.0) RETURN XC CLEAR THE ERROR TABLES X DO 70 I=1,10 X 70 KOUNT(I) = 0 X KOUNTX = 0 X RETURN X 80 CONTINUE XC PROCESS A MESSAGE... XC SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG, XC OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL. X MES = MESSG X DO 90 I=1,10 X II = I X IF (KOUNT(I).EQ.0) GO TO 110 X IF (MES.NE.MESTAB(I)) GO TO 90 X IF (NERR.NE.NERTAB(I)) GO TO 90 X IF (LEVEL.NE.LEVTAB(I)) GO TO 90 X GO TO 100 X 90 CONTINUE XC THREE POSSIBLE CASES... XC TABLE IS FULL X KOUNTX = KOUNTX+1 X ICOUNT = 1 X RETURN XC MESSAGE FOUND IN TABLE X 100 KOUNT(II) = KOUNT(II) + 1 X ICOUNT = KOUNT(II) X RETURN XC EMPTY SLOT FOUND FOR NEW MESSAGE X 110 MESTAB(II) = MES X NERTAB(II) = NERR X LEVTAB(II) = LEVEL X KOUNT(II) = 1 X ICOUNT = 1 X RETURN X END X SUBROUTINE XGETF(KONTRL) XC***BEGIN PROLOGUE XGETF XC***DATE WRITTEN 790801 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3C XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Returns current value of error control flag. XC***DESCRIPTION XC Abstract XC XGETF returns the current value of the error control flag XC in KONTRL. See subroutine XSETF for flag value meanings. XC (KONTRL is an output parameter only.) XC XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC Latest revision --- 7 June 1978 XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED J4SAVE XC***END PROLOGUE XGETF XC***FIRST EXECUTABLE STATEMENT XGETF X INTEGER J4SAVE X INTEGER KONTRL X KONTRL = J4SAVE(2,0,.FALSE.) X RETURN X END X SUBROUTINE XGETUA(IUNITA,N) XC***BEGIN PROLOGUE XGETUA XC***DATE WRITTEN 790801 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3C XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Returns unit number(s) to which error messages are being XC sent. XC***DESCRIPTION XC Abstract XC XGETUA may be called to determine the unit number or numbers XC to which error messages are being sent. XC These unit numbers may have been set by a call to XSETUN, XC or a call to XSETUA, or may be a default value. XC XC Description of Parameters XC --Output-- XC IUNIT - an array of one to five unit numbers, depending XC on the value of N. A value of zero refers to the XC default unit, as defined by the I1MACH machine XC constant routine. Only IUNIT(1),...,IUNIT(N) are XC defined by XGETUA. The values of IUNIT(N+1),..., XC IUNIT(5) are not defined (for N .LT. 5) or altered XC in any way by XGETUA. XC N - the number of units to which copies of the XC error messages are being sent. N will be in the XC range from 1 to 5. XC XC Latest revision --- 19 MAR 1980 XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED J4SAVE XC***END PROLOGUE XGETUA X INTEGER J4SAVE X INTEGER I, INDEX, IUNITA, N X DIMENSION IUNITA(5) XC***FIRST EXECUTABLE STATEMENT XGETUA X N = J4SAVE(5,0,.FALSE.) X DO 30 I=1,N X INDEX = I+4 X IF (I.EQ.1) INDEX = 3 X IUNITA(I) = J4SAVE(INDEX,0,.FALSE.) X 30 CONTINUE X RETURN X END X SUBROUTINE XSETF(KONTRL) XC***BEGIN PROLOGUE XSETF XC***DATE WRITTEN 790801 (YYMMDD) XC***REVISION DATE 820801 (YYMMDD) XC***CATEGORY NO. R3A XC***KEYWORDS ERROR,XERROR PACKAGE XC***AUTHOR JONES, R. E., (SNLA) XC***PURPOSE Sets the error control flag. XC***DESCRIPTION XC Abstract XC XSETF sets the error control flag value to KONTRL. XC (KONTRL is an input parameter only.) XC The following table shows how each message is treated, XC depending on the values of KONTRL and LEVEL. (See XERROR XC for description of LEVEL.) XC XC If KONTRL is zero or negative, no information other than the XC message itself (including numeric values, if any) will be XC printed. If KONTRL is positive, introductory messages, XC trace-backs, etc., will be printed in addition to the message. XC XC IABS(KONTRL) XC LEVEL 0 1 2 XC value XC 2 fatal fatal fatal XC XC 1 not printed printed fatal XC XC 0 not printed printed printed XC XC -1 not printed printed printed XC only only XC once once XC XC Written by Ron Jones, with SLATEC Common Math Library Subcommittee XC Latest revision --- 19 MAR 1980 XC***REFERENCES JONES R.E., KAHANER D.K., "XERROR, THE SLATEC ERROR- XC HANDLING PACKAGE", SAND82-0800, SANDIA LABORATORIES, XC 1982. XC***ROUTINES CALLED J4SAVE,XERRWV XC***END PROLOGUE XSETF XC***FIRST EXECUTABLE STATEMENT XSETF X INTEGER J4SAVE X INTEGER JUNK, KONTRL X IF ((KONTRL.GE.(-2)).AND.(KONTRL.LE.2)) GO TO 10 X CALL XERRWV('XSETF -- INVALID VALUE OF KONTRL (I1).',33,1,2, X 1 1,KONTRL,0,0,0.,0.) X RETURN X 10 JUNK = J4SAVE(2,KONTRL,.TRUE.) X RETURN X END X DOUBLE PRECISION FUNCTION D1MACH(I) XC***BEGIN PROLOGUE D1MACH XC***DATE WRITTEN 750101 (YYMMDD) XC***REVISION DATE 910131 (YYMMDD) XC***CATEGORY NO. R1 XC***KEYWORDS MACHINE CONSTANTS XC***AUTHOR FOX, P. A., (BELL LABS) XC HALL, A. D., (BELL LABS) XC SCHRYER, N. L., (BELL LABS) XC***PURPOSE Returns double precision machine dependent constants XC***DESCRIPTION XC XC This is the CMLIB version of D1MACH, the double precision machine XC constants subroutine originally developed for the PORT library. XC XC D1MACH can be used to obtain machine-dependent parameters XC for the local machine environment. It is a function XC subprogram with one (input) argument, and can be called XC as follows, for example XC XC D = D1MACH(I) XC XC where I=1,...,5. The (output) value of D above is XC determined by the (input) value of I. The results for XC various values of I are discussed below. XC XC Double-precision machine constants XC D1MACH( 1) = B**(EMIN-1), the smallest positive magnitude. XC D1MACH( 2) = B**EMAX*(1 - B**(-T)), the largest magnitude. XC D1MACH( 3) = B**(-T), the smallest relative spacing. XC D1MACH( 4) = B**(1-T), the largest relative spacing. XC D1MACH( 5) = LOG10(B) XC***REFERENCES FOX P.A., HALL A.D., SCHRYER N.L.,*FRAMEWORK FOR A XC PORTABLE LIBRARY*, ACM TRANSACTIONS ON MATHEMATICAL XC SOFTWARE, VOL. 4, NO. 2, JUNE 1978, PP. 177-188. XC***ROUTINES CALLED XERROR XC***END PROLOGUE D1MACH XC X INTEGER I X INTEGER SMALL(4) X INTEGER LARGE(4) X INTEGER RIGHT(4) X INTEGER DIVER(4) X INTEGER LOG10(4) XC X DOUBLE PRECISION DMACH(5) XC X EQUIVALENCE (DMACH(1),SMALL(1)) X EQUIVALENCE (DMACH(2),LARGE(1)) X EQUIVALENCE (DMACH(3),RIGHT(1)) X EQUIVALENCE (DMACH(4),DIVER(1)) X EQUIVALENCE (DMACH(5),LOG10(1)) XC XC MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES, SUCH AS THE AT&T XC 3B SERIES AND MOTOROLA 68000 BASED MACHINES (E.G. SUN 3 AND AT&T XC PC 7300), IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST. XC XC === MACHINE = IEEE.MOST-SIG-BYTE-FIRST XC === MACHINE = SUN XC === MACHINE = 68000 XC === MACHINE = ATT.3B XC === MACHINE = ATT.7300 X DATA SMALL(1),SMALL(2) / 1048576, 0 / X DATA LARGE(1),LARGE(2) / 2146435071, -1 / X DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / X DATA DIVER(1),DIVER(2) / 1018167296, 0 / X DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 / XC XC MACHINE CONSTANTS FOR IEEE ARITHMETIC MACHINES AND 8087-BASED XC MICROS, SUCH AS THE IBM PC AND AT&T 6300, IN WHICH THE LEAST XC SIGNIFICANT BYTE IS STORED FIRST. XC XC === MACHINE = IEEE.LEAST-SIG-BYTE-FIRST XC === MACHINE = 8087 XC === MACHINE = IBM.PC XC === MACHINE = ATT.6300 XC DATA SMALL(1),SMALL(2) / 0, 1048576 / XC DATA LARGE(1),LARGE(2) / -1, 2146435071 / XC DATA RIGHT(1),RIGHT(2) / 0, 1017118720 / XC DATA DIVER(1),DIVER(2) / 0, 1018167296 / XC DATA LOG10(1),LOG10(2) / 1352628735, 1070810131 / XC XC MACHINE CONSTANTS FOR AMDAHL MACHINES. XC XC === MACHINE = AMDAHL XC DATA SMALL(1),SMALL(2) / 1048576, 0 / XC DATA LARGE(1),LARGE(2) / 2147483647, -1 / XC DATA RIGHT(1),RIGHT(2) / 856686592, 0 / XC DATA DIVER(1),DIVER(2) / 873463808, 0 / XC DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 / XC XC MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. XC XC === MACHINE = BURROUGHS.1700 XC DATA SMALL(1) / ZC00800000 / XC DATA SMALL(2) / Z000000000 / XC DATA LARGE(1) / ZDFFFFFFFF / XC DATA LARGE(2) / ZFFFFFFFFF / XC DATA RIGHT(1) / ZCC5800000 / XC DATA RIGHT(2) / Z000000000 / XC DATA DIVER(1) / ZCC6800000 / XC DATA DIVER(2) / Z000000000 / XC DATA LOG10(1) / ZD00E730E7 / XC DATA LOG10(2) / ZC77800DC0 / XC XC MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. XC XC === MACHINE = BURROUGHS.5700 XC DATA SMALL(1) / O1771000000000000 / XC DATA SMALL(2) / O0000000000000000 / XC DATA LARGE(1) / O0777777777777777 / XC DATA LARGE(2) / O0007777777777777 / XC DATA RIGHT(1) / O1461000000000000 / XC DATA RIGHT(2) / O0000000000000000 / XC DATA DIVER(1) / O1451000000000000 / XC DATA DIVER(2) / O0000000000000000 / XC DATA LOG10(1) / O1157163034761674 / XC DATA LOG10(2) / O0006677466732724 / XC XC MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. XC XC === MACHINE = BURROUGHS.6700 XC === MACHINE = BURROUGHS.7700 XC DATA SMALL(1) / O1771000000000000 / XC DATA SMALL(2) / O7770000000000000 / XC DATA LARGE(1) / O0777777777777777 / XC DATA LARGE(2) / O7777777777777777 / XC DATA RIGHT(1) / O1461000000000000 / XC DATA RIGHT(2) / O0000000000000000 / XC DATA DIVER(1) / O1451000000000000 / XC DATA DIVER(2) / O0000000000000000 / XC DATA LOG10(1) / O1157163034761674 / XC DATA LOG10(2) / O0006677466732724 / XC XC MACHINE CONSTANTS FOR THE CONVEX C-120 (NATIVE MODE) XC WITH OR WITHOUT -R8 OPTION XC XC === MACHINE = CONVEX.C1 XC === MACHINE = CONVEX.C1.R8 XC DATA DMACH(1) / 5.562684646268007D-309 / XC DATA DMACH(2) / 8.988465674311577D+307 / XC DATA DMACH(3) / 1.110223024625157D-016 / XC DATA DMACH(4) / 2.220446049250313D-016 / XC DATA DMACH(5) / 3.010299956639812D-001 / XC XC MACHINE CONSTANTS FOR THE CONVEX C-120 (IEEE MODE) XC WITH OR WITHOUT -R8 OPTION XC XC === MACHINE = CONVEX.C1.IEEE XC === MACHINE = CONVEX.C1.IEEE.R8 XC DATA DMACH(1) / 2.225073858507202D-308 / XC DATA DMACH(2) / 1.797693134862315D+308 / XC DATA DMACH(3) / 1.110223024625157D-016 / XC DATA DMACH(4) / 2.220446049250313D-016 / XC DATA DMACH(5) / 3.010299956639812D-001 / XC XC MACHINE CONSTANTS FOR THE CYBER 170/180 SERIES USING NOS (FTN5). XC XC === MACHINE = CYBER.170.NOS XC === MACHINE = CYBER.180.NOS XC DATA SMALL(1) / O"00604000000000000000" / XC DATA SMALL(2) / O"00000000000000000000" / XC DATA LARGE(1) / O"37767777777777777777" / XC DATA LARGE(2) / O"371677