SUBROUTINE VJ1 (M, X, F, WORK, IWORK, INFO) C***BEGIN PROLOGUE VJ1 C***PURPOSE Computes the Bessel function of the first kind C of order one (J1) for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C10A1 C***TYPE SINGLE PRECISION (VJ1-S, DVJ1-D) C***KEYWORDS BESSEL FUNCTION, FIRST KIND, SPECIAL FUNCTION, C ORDER ONE, VECTORIZED C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C VJ1 computes the Bessel function of the first kind of order C one (J1) for real arguments using uniform approximation by C Chebyshev polynomials. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of arguments at which the function is to be C evaluated. C C X (Input) Real array of length M C The arguments at which the function is to be evaluated are C stored in X(1) to X(M) in any order. C C F (Output) Real array of length M C F(i) contains the value of the function at X(i), i=1,..,M. C C WORK (Work) Real vector of length 7*M C C IWORK (Work) Integer vector of length M C C C INFO (Output) Integer C Indicates status of computed result. The following table C lists possible values and their meanings. If OK=Yes then C all F(i) have been set, otherwise none have been set. C C INFO OK Description C ------------------------------------------------------------ C -1 Yes Warning: Some abs(X(i)) so small J1 underflows. C The corresponding F(i) are set to zero. C 0 Yes Successfull execution. C 1 No Error: M .LE. 0 C 2 No Error: Some abs(X(i)) so big that no precision C possible in computing J1. The index of the C first offending argument is returned in IWORK(1). C C ********************************************************************* C This routine is a modification of the function BESJ1 developed by C W. Fullerton of LANL. C ********************************************************************* C C***REFERENCES Ronald F. Boisvert and Bonita V. Saunders, Portable C Vectorized Software for Bessel Function Evaluation, C ACM Transactions on Mathematical Software 18 (1992), C pp 456-469. C***ROUTINES CALLED WJ1 C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE VJ1 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, IWORK, M REAL F, X, WORK C DIMENSION X(M), F(M), WORK(7*M), IWORK(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER IWB0, IWB1, IWB2, IWY, IWTC, IWYC, IWZC, JIN C C***FIRST EXECUTABLE STATEMENT VJ1 C C ... PARTITION WORK ARRAYS C IWY = 1 IWTC = IWY + M IWYC = IWTC + M IWZC = IWYC + M IWB0 = IWZC + M IWB1 = IWB0 + M IWB2 = IWB1 + M C Total = IWB2 + M C JIN = 1 C Total = JIN + M C C ... WJ1 DOES ALL THE WORK C CALL WJ1(M,X,F,WORK(IWY),WORK(IWTC),WORK(IWYC),WORK(IWZC), + IWORK(JIN),WORK(IWB0),WORK(IWB1),WORK(IWB2),INFO) C RETURN END SUBROUTINE WJ1(M, X, F, Y, TCMP, YCMP, ZCMP, INDX, B0, B1, B2, + INFO) C***BEGIN PROLOGUE WJ1 C***SUBSIDIARY C***PURPOSE Computes the Bessel function of the first kind C of order one (J1) for a vector of arguments C***LIBRARY VFNLIB C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***ROUTINES CALLED IWCS, R1MACH, WNGT, WGTHR, WGTLE, WGT, C WLE, WSCTR, WCS C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WJ1 C C ---------- C PARAMETERS C ---------- C INTEGER INFO, INDX, M REAL B0, B1, B2, F, X, Y, TCMP, YCMP, ZCMP C DIMENSION B0(M), B1(M), B2(M), F(M), INDX(M), X(M), Y(M), + TCMP(M), YCMP(M), ZCMP(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER LBJ1, LBM1, LBM12, LBTH1, LBT12 PARAMETER (LBJ1=12, LBM1=21, LBM12=2, LBTH1=24, LBT12=2) C INTEGER I, IWCS, J, JH, K, KEY, N, NA, NB, NTJ1, NTM1, NTM12, + NTTH1, NTT12, NTOT REAL BJ1CS, BM1CS, BM12CS, BTH1CS, BT12CS, C1, C2, + EPMACH, EPS, PI4, R1MACH, TPI4, XMAX, XMIN, XSML C DIMENSION BJ1CS(LBJ1), BM1CS(LBM1), BM12CS(LBM12), BTH1CS(LBTH1), + BT12CS(LBT12) C SAVE BJ1CS, BM1CS, BM12CS, BTH1CS, BT12CS, NTJ1, NTM1, NTM12, + NTTH1, NTT12, PI4, TPI4, XMIN, XSML, XMAX C C---------------------------------------------------------------------- C C Series for BJ1 on the interval 0. to 1.60000D+01 C with weighted error 4.48E-17 C log weighted error 16.35 C significant figures required 15.77 C decimal places required 16.89 C DATA BJ1 CS( 1) / -.1172614151 3332787E0 / DATA BJ1 CS( 2) / -.2536152183 0790640E0 / DATA BJ1 CS( 3) / .0501270809 84469569E0 / DATA BJ1 CS( 4) / -.0046315148 09625081E0 / DATA BJ1 CS( 5) / .0002479962 29415914E0 / DATA BJ1 CS( 6) / -.0000086789 48686278E0 / DATA BJ1 CS( 7) / .0000002142 93917143E0 / DATA BJ1 CS( 8) / -.0000000039 36093079E0 / DATA BJ1 CS( 9) / .0000000000 55911823E0 / DATA BJ1 CS(10) / -.0000000000 00632761E0 / DATA BJ1 CS(11) / .0000000000 00005840E0 / DATA BJ1 CS(12) / -.0000000000 00000044E0 / C C---------------------------------------------------------------------- C C Series for BM1 on the interval 0. to 6.25000D-02 C with weighted error 5.61E-17 C log weighted error 16.25 C significant figures required 14.97 C decimal places required 16.91 C DATA BM1 CS( 1) / .1047362510 931285E0 / DATA BM1 CS( 2) / .0044244389 3702345E0 / DATA BM1 CS( 3) / -.0000566163 9504035E0 / DATA BM1 CS( 4) / .0000023134 9417339E0 / DATA BM1 CS( 5) / -.0000001737 7182007E0 / DATA BM1 CS( 6) / .0000000189 3209930E0 / DATA BM1 CS( 7) / -.0000000026 5416023E0 / DATA BM1 CS( 8) / .0000000004 4740209E0 / DATA BM1 CS( 9) / -.0000000000 8691795E0 / DATA BM1 CS(10) / .0000000000 1891492E0 / DATA BM1 CS(11) / -.0000000000 0451884E0 / DATA BM1 CS(12) / .0000000000 0116765E0 / DATA BM1 CS(13) / -.0000000000 0032265E0 / DATA BM1 CS(14) / .0000000000 0009450E0 / DATA BM1 CS(15) / -.0000000000 0002913E0 / DATA BM1 CS(16) / .0000000000 0000939E0 / DATA BM1 CS(17) / -.0000000000 0000315E0 / DATA BM1 CS(18) / .0000000000 0000109E0 / DATA BM1 CS(19) / -.0000000000 0000039E0 / DATA BM1 CS(20) / .0000000000 0000014E0 / DATA BM1 CS(21) / -.0000000000 0000005E0 / C C---------------------------------------------------------------------- C C Series for BM12 Used in double precision version only C DATA BM12CS( 1) / 1.0E0 / DATA BM12CS( 2) / 0.0E0 / C C---------------------------------------------------------------------- C C Series for BTH1 on the interval 0. to 6.25000D-02 C with weighted error 4.10E-17 C log weighted error 16.39 C significant figures required 15.96 C DATA BTH1CS( 1) / .7406014102 6313850E0 / DATA BTH1CS( 2) / -.0045717556 59637690E0 / DATA BTH1CS( 3) / .0001198185 10964326E0 / DATA BTH1CS( 4) / -.0000069645 61891648E0 / DATA BTH1CS( 5) / .0000006554 95621447E0 / DATA BTH1CS( 6) / -.0000000840 66228945E0 / DATA BTH1CS( 7) / .0000000133 76886564E0 / DATA BTH1CS( 8) / -.0000000024 99565654E0 / DATA BTH1CS( 9) / .0000000005 29495100E0 / DATA BTH1CS(10) / -.0000000001 24135944E0 / DATA BTH1CS(11) / .0000000000 31656485E0 / DATA BTH1CS(12) / -.0000000000 08668640E0 / DATA BTH1CS(13) / .0000000000 02523758E0 / DATA BTH1CS(14) / -.0000000000 00775085E0 / DATA BTH1CS(15) / .0000000000 00249527E0 / DATA BTH1CS(16) / -.0000000000 00083773E0 / DATA BTH1CS(17) / .0000000000 00029205E0 / DATA BTH1CS(18) / -.0000000000 00010534E0 / DATA BTH1CS(19) / .0000000000 00003919E0 / DATA BTH1CS(20) / -.0000000000 00001500E0 / DATA BTH1CS(21) / .0000000000 00000589E0 / DATA BTH1CS(22) / -.0000000000 00000237E0 / DATA BTH1CS(23) / .0000000000 00000097E0 / DATA BTH1CS(24) / -.0000000000 00000040E0 / C C---------------------------------------------------------------------- C C Series for BT12 Used in double precision version only C DATA BT12CS( 1) / 1.0E0 / DATA BT12CS( 2) / 0.0E0 / C C---------------------------------------------------------------------- C DATA NTJ1 / 0 / C C***FIRST EXECUTABLE STATEMENT WJ1 C IF (M .LE. 0) GO TO 910 C IF (NTJ1.EQ.0) THEN EPMACH = R1MACH(3) EPS = 0.10E0*EPMACH NTJ1 = IWCS(BJ1CS, LBJ1, EPS) NTM1 = IWCS(BM1CS, LBM1, EPS) NTM12 = IWCS(BM12CS, LBM12, EPS) NTTH1 = IWCS(BTH1CS, LBTH1, EPS) NTT12 = IWCS(BT12CS, LBT12, EPS) XMIN = 2.0E0*R1MACH(1) XSML = SQRT(8.0E0*EPMACH) XMAX = 1.0E0/R1MACH(4) PI4 = ATAN(1.0E0) TPI4 = 3.0E0*PI4 ENDIF C NTOT = 0 C DO 10 I=1,M Y(I) = ABS(X(I)) 10 CONTINUE C CALL WNGT(M,Y,XMAX,KEY) IF (KEY .NE. 0) GO TO 920 C C --------------------------- C CASE Y=0 OR Y TOO SMALL C --------------------------- C C NOTE --- J1 UNDERFLOWS FOR 0 .LT. Y .LE. XMIN C DO 15 I=1,M F(I) = 0.0E0 15 CONTINUE C C ---------------------------- C CASE XMIN .LT. Y .LE. XSML C ---------------------------- C CALL WGTLE(M,Y,XMIN,XSML,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 20 J=1,N ZCMP(J) = 0.50E0*YCMP(J) 20 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C --------------------------- C CASE XSML .LT. Y .LE. 4.0 C --------------------------- C CALL WGTLE(M,Y,XSML,4.0E0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 30 J=1,N TCMP(J) = 0.1250E0*YCMP(J)**2 - 1.0E0 30 CONTINUE CALL WCS(N,TCMP,BJ1CS,NTJ1,ZCMP,B0,B1,B2) DO 40 K=1,N ZCMP(K) = YCMP(K)*(0.250E0 + ZCMP(K)) 40 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ---------------- C CASE Y .GT. 4.0 C ---------------- C CALL WGT(M,Y,4.0E0,N,INDX) IF (N .GT. 0) THEN NTOT = NTOT + N CALL WGTHR(N,Y,INDX,YCMP) DO 50 J=1,N TCMP(J) = 32.0E0/YCMP(J)**2 - 1.0E0 50 CONTINUE CALL WCS(N,TCMP,BM1CS,NTM1,Y,B0,B1,B2) CALL WCS(N,TCMP,BTH1CS,NTTH1,ZCMP,B0,B1,B2) DO 60 J=1,N Y(J) = (0.750E0 + Y(J)) / SQRT(YCMP(J)) ZCMP(J) = ( YCMP(J) - TPI4 ) + ZCMP(J)/YCMP(J) ZCMP(J) = Y(J) * COS(ZCMP(J)) 60 CONTINUE CALL WSCTR(N,ZCMP,INDX,F) ENDIF C C ------------------------------------------- C REVERSE SIGN FOR NEGATIVE X (RESULT IS ODD) C ------------------------------------------- C DO 200 I=1,M IF (X(I) .LT. 0.0E0) F(I) = -F(I) 200 CONTINUE C C ----- C EXITS C ----- C C ... NORMAL C INFO = 0 IF (NTOT .NE. M) INFO = -1 GO TO 999 C C ... M .LE. 0 C 910 CONTINUE INFO = 1 GO TO 999 C C ... NO PRECISION BECAUSE X BIG C 920 CONTINUE INFO = 2 INDX(1) = KEY GO TO 999 C 999 CONTINUE RETURN END SUBROUTINE WCS (M, X, CS, N, F, B0, B1, B2) C***BEGIN PROLOGUE WCS C***PURPOSE Evaluate a Chebyshev series for a vector of real arguments C***LIBRARY VFNLIB C***CATEGORY C3A2 C***TYPE SINGLE PRECISION (WCS-S, DWCS-D) C***KEYWORDS CHEBYSHEV SERIES EVALUATION, VECTORIZED C***AUTHOR SAUNDERS, B. V,. (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WCS evaluates a given Chebyshev series for a vector of real C arguments. C C C P A R A M E T E R S C C M (Input) Integer C The number of arguments at which the series is to be C evaluated. C C X (Input) Real array of length M C The arguments at which the function is to be evaluated are C stored in X(1) to X(M) in any order. C C CS (Input) Real array of length .GE. N C The N coefficients of the Chebyshev series. (Note -- only half C the first coefficient is used in summing the series.) C C N (Input) Integer (0 .LE. N .LE. 1000) C The number of terms in the Chebyshev series. C C F (Output) Real array of length M C F(i) contains the value of the series at X(i), i=1,..,M. C C B0 (Work) Real array of length M C C B1 (Work) Real array of length M C C B2 (Work) Real array of length M C C C ********************************************************************* C This routine is a modification of the function CSEVL developed by C W. Fullerton of LANL. C ********************************************************************* C C***REFERENCES Ronald F. Boisvert and Bonita V. Saunders, Portable C Vectorized Software for Bessel Function Evaluation, C ACM Transactions on Mathematical Software 18 (1992), C pp 456-469. C***ROUTINES CALLED WFERR C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE WCS C C ---------- C PARAMETERS C ---------- C INTEGER N, M REAL X, CS, F, B0, B1, B2 C DIMENSION B0(M), B1(M), B2(M), CS(N), F(M), X(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I, IMOD, J, NI REAL CS1, CS2, CSNI, CSNI1, CSNI2 C C***FIRST EXECUTABLE STATEMENT WCS C IF (N .LT. 0) CALL WFERR('WCS','NUMBER OF TERMS LT 0',2) IF (N .GT. 1000) CALL WFERR('WCS','NUMBER OF TERMS GT 1000',2) C C ... INITIALIZATION C DO 10 I=1,M B1(I) = 0.0E0 B2(I) = 0.0E0 F(I) = 2.0E0*X(I) 10 CONTINUE C C ... THREE-TERM RECURSION C (DO THREE STEPS AT A TIME TO AVOID VECTOR COPIES) C IMOD = MOD(N,3) DO 30 I=1,N-IMOD,3 NI = N + 1 - I CSNI = CS(NI) CSNI1 = CS(NI-1) CSNI2 = CS(NI-2) DO 20 J=1,M B0(J) = ( F(J)*B1(J)-B2(J) ) + CSNI B2(J) = ( F(J)*B0(J)-B1(J) ) + CSNI1 B1(J) = ( F(J)*B2(J)-B0(J) ) + CSNI2 20 CONTINUE 30 CONTINUE C C ... LAST STEP C (CLEANUP FOR CASE N NOT DIVISIBLE BY 3) C IF (IMOD .EQ. 0) THEN DO 40 J=1,M F(J) = 0.50E0*(B1(J) - B0(J)) 40 CONTINUE ELSEIF (IMOD .EQ. 1) THEN CS1 = CS(1) DO 50 J=1,M B0(J) = ( F(J)*B1(J) - B2(J) ) + CS1 F(J) = 0.50E0*(B0(J) - B2(J)) 50 CONTINUE ELSE CS1 = CS(1) CS2 = CS(2) DO 60 J=1,M B0(J) = ( F(J)*B1(J) - B2(J) ) + CS2 B2(J) = ( F(J)*B0(J) - B1(J) ) + CS1 F(J) = 0.50E0*(B2(J) - B1(J)) 60 CONTINUE ENDIF C RETURN END INTEGER FUNCTION IWCS (OS, NOS, ETA) C***BEGIN PROLOGUE IWCS C***PURPOSE Determines the number of terms of an orthogonal series C required to meet a specified error tolerance. C***LIBRARY VFNLIB C***CATEGORY C3A2 C***TYPE SINGLE PRECISION (IWCS-S, IDWCS-D) C***KEYWORDS CHEBYSHEV SERIES, INITIALIZATION C***AUTHOR BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C IWCS returns the number of terms of the Chebyshev series OS C required to insure the error in evaluating it is no larger than C ETA. C C Ordinarily, ETA is chosen to be one-tenth machine precision. C C C P A R A M E T E R S C C OS (Input) Real array of length .GE. NOS C Coefficients of the NOS-term orthogonal series. C C NOS (Input) Integer (.GE. 1) C Number of terms in the orthogonal series. C C ETA (Input) Real C Requested accuracy of the series. C C C ********************************************************************* C This routine is a modification of the function INITS developed by C W. Fullerton of LANL. C ********************************************************************* C C***REFERENCES Ronald F. Boisvert and Bonita V. Saunders, Portable C Vectorized Software for Bessel Function Evaluation, C ACM Transactions on Mathematical Software 18 (1992), C pp 456-469. C***ROUTINES CALLED WFERR C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C 930203 Minor modifications to prologue. C***END PROLOGUE IWCS C C ---------- C PARAMETERS C ---------- C INTEGER NOS REAL ETA, OS C DIMENSION OS(NOS) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I, II REAL ERR C C***FIRST EXECUTABLE STATEMENT IWCS C IF (NOS .LT. 1) + CALL WFERR('IWCS','NUMBER OF COEFFICIENTS LT 1',2) C ERR = 0.0E0 DO 10 II=NOS,1,-1 I = II ERR = ERR + ABS(OS(I)) IF (ERR .GT. ETA) GO TO 20 10 CONTINUE C 20 CONTINUE IF (I .EQ. NOS) + CALL WFERR('IWCS','TOO MUCH ACCURACY REQUESTED',2) C IWCS = I C RETURN END SUBROUTINE WFERR (SNAME, MESSAG, LEVEL) C***BEGIN PROLOGUE WFERR C***PURPOSE Processes an error message for the vectorized special C function package C***LIBRARY VFNLIB C***TYPE CHARACTER C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WFERR processes an error message for the vectorized special C function package. C C P A R A M E T E R S C C SNAME (Input) Character*(*) C The name of the routine issuing the message. C C MESSAG (Input) Character*(*) C The error message. Must be less than 80 characters. C C LEVEL (Input) Integer C The message level. 1 for recoverable, 2 for fatal. C C This routine issues a STOP for fatal errors. C C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WFERR C C ... PARAMETERS C CHARACTER*(*) MESSAG, SNAME INTEGER LEVEL C C ... LOCAL VARIABLES C INTEGER IUNIT, I1MACH, MLEN C C***FIRST EXECUTABLE STATEMENT WFERR C IUNIT = I1MACH(4) MLEN = MIN(LEN(MESSAG),80) C WRITE(IUNIT,*) IF (LEVEL .EQ. 2) THEN WRITE(IUNIT,*) 'FATAL ERROR IN ',SNAME,' ...' WRITE(IUNIT,*) MESSAG(1:MLEN) STOP ELSE WRITE(IUNIT,*) 'RECOVERABLE ERROR IN ',SNAME,' ...' WRITE(IUNIT,*) MESSAG(1:MLEN) RETURN ENDIF C END SUBROUTINE WGT (M, X, SCALR, N, INDX) C***BEGIN PROLOGUE WGT C***PURPOSE Builds an array of indices of vector elements that are C greater than a specified scalar C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WGT-S, DWGT-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WGT builds an array of indices of vector elements that are C greater than a specified scalar. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of elements of the input vector. C C X (Input) Real array of length M C The input vector. C C SCALR (Input) Real C The scalar used for comparison with vector elements. C Indices i are selected if X(i) .GT. SCALR. C C N (Output) Integer C The number of elements of X that satisfy relationship with C scalar. C C INDX (Output) Integer array of length N C Array containing indices of vector elements that satisfy C relationship with scalar. C C C CAUTION : constraints on the input variables are not checked by C this routine. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WGT C C ---------- C PARAMETERS C ---------- C INTEGER INDX, M, N REAL SCALR, X C DIMENSION X(M), INDX(*) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I C C***FIRST EXECUTABLE STATEMENT WGT C N = 0 C DO 10 I=1,M IF (X(I) .GT. SCALR) THEN N = N + 1 INDX(N) = I ENDIF 10 CONTINUE C RETURN END SUBROUTINE WGTHR (N, X, INDX, Y) C***BEGIN PROLOGUE WGTHR C***PURPOSE Performs a gather by index on a given vector C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WGTHR-S, DWGTHR-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WGTHR performs a gather by index on a given vector based on the C indices provided in an array: Y(i) = X(INDX(i)), i=1,...,N. C C C P A R A M E T E R S C C N (Input) Integer (N .GT. 0) C The number of elements to be gathered from the input vector C X. C C X (Input) Real array of length M C The input vector. C C INDX (Input) Integer array of length N C Array specifying indices of elements to be gathered from C the input vector. C C Y (Output) Real array of length N C The array containing the compressed vector once the gather is C completed. Y(i) = X(INDX(i)), i=1,...,N. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WGTHR C C ---------- C PARAMETERS C ---------- C INTEGER INDX, N REAL X, Y C DIMENSION X(*), Y(N), INDX(N) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER J C C***FIRST EXECUTABLE STATEMENT WGTHR C IF (N .LE. 0) RETURN C DO 10 J=1,N Y(J) = X(INDX(J)) 10 CONTINUE C RETURN END SUBROUTINE WGTLE (M, X, SCALR1, SCALR2, N, INDX) C***BEGIN PROLOGUE WGTLE C***PURPOSE Builds an array of indices of vector elements that are C greater than a given scalar and less than or equal to a C second given scalar C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WGTLE-S, DWGTLE-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WGTLE builds an array of indices of vector elements that are C greater than a given scalar and less than or equal to a second C given scalar. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of elements of the input vector. C C X (Input) Real array of length M C The input vector. C C SCALR1, C SCALR2 (Input) Real C The scalars used for comparison with vector elements. C Indices i are selected if SCALR1 .LT. X(i) .LE. SCALR2. C C N (Output) Integer C The number of elements of X that satisfy relationship with C scalars. C C INDX (Output) Integer array of length N C Array containing indices of vector elements that satisfy C relationship with scalars. C C C CAUTION : constraints on the input variables are not checked by C this routine. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WGTLE C C ---------- C PARAMETERS C ---------- C INTEGER INDX, M, N REAL SCALR1, SCALR2, X C DIMENSION X(M), INDX(*) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I REAL ELEMT C C***FIRST EXECUTABLE STATEMENT WGTLE C N = 0 C DO 10 I=1,M ELEMT = X(I) IF (ELEMT .GT. SCALR1 .AND. ELEMT .LE. SCALR2) THEN N = N + 1 INDX(N) = I ENDIF 10 CONTINUE C RETURN END SUBROUTINE WLE (M, X, SCALR, N, INDX) C***BEGIN PROLOGUE WLE C***PURPOSE Builds an array of indices of vector elements that are C less than or equal to a given scalar C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WLE-S, DWLE-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WLE builds an array of indices of vector elements that are C less than or equal to a given scalar. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of elements of the input vector. C C X (Input) Real array of length M C The input vector. C C SCALR (Input) Real C The scalar used for comparison with vector elements. C Indices i are selected if X(i) .LE. SCALR. C C N (Output) Integer C The number of elements of X that satisfy relationship with C scalars. C C INDX (Output) Integer array of length N C Array containing indices of vector elements that satisfy C relationship with scalars. C C C CAUTION : constraints on the input variables are not checked by C this routine. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WLE C C ---------- C PARAMETERS C ---------- C INTEGER INDX, M, N REAL SCALR, X C DIMENSION X(M), INDX(*) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I C C***FIRST EXECUTABLE STATEMENT WGE C N = 0 C DO 10 I = 1,M IF (X(I) .LE. SCALR) THEN N = N + 1 INDX(N) = I ENDIF 10 CONTINUE C RETURN END SUBROUTINE WNGT (M, X, SCALR, KEY) C***BEGIN PROLOGUE WNGT C***PURPOSE Determines whether elements of a given vector are greater C than a given scalar C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WNGT-S, DWNGT-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WNGT determines whether elements of a given vector are greater C than a given scalar. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of elements of the input vector. C C X (Input) Real array of length M C The input vector. C C SCALR (Input) Real C The scalar used for comparison with vector elements. C C KEY (Output) Integer C If KEY=0 then no vector elements satisfy X(i).GT.SCALR. C If KEY>0 then some vector elements satisfy X(i).GT.SCALR; C the first to do so is X(KEY). C C C CAUTION : constraints on the input variables are not checked by C this routine. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WNGT C C ---------- C PARAMETERS C ---------- C INTEGER KEY, M REAL SCALR, X C DIMENSION X(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I C C***FIRST EXECUTABLE STATEMENT WNGT C C ... QUICK CHECK (VECTORIZABLE) C KEY = 0 DO 10 I=1,M IF (X(I) .GT. SCALR) KEY = 1 10 CONTINUE C C ... IF CHECK FAILED THEN FIND INDEX OF FIRST FAILURE C IF (KEY .NE. 0) THEN KEY = 0 DO 20 I=1,M IF (X(I) .GT. SCALR) THEN KEY = I GO TO 30 ENDIF 20 CONTINUE ENDIF C 30 CONTINUE RETURN END SUBROUTINE WNLE (M, X, SCALR, KEY) C***BEGIN PROLOGUE WNLE C***PURPOSE Determines whether elements of a given vector are less C than or equal to a given scalar C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WNLE-S, DWNLE-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WNLE determines whether elements of a given vector are less than C or equal to a given scalar. C C C P A R A M E T E R S C C M (Input) Integer (M .GT. 0) C The number of elements of the input vector. C C X (Input) Real array of length M C The input vector. C C SCALR (Input) Real C The scalar used for comparison with vector elements. C C KEY (Output) Integer C If KEY=0 then no vector elements satisfy X(i).LE.SCALR. C If KEY>0 then some vector elements satisfy X(i).LE.SCALR; C the first to do so is X(KEY). C C C CAUTION : constraints on the input variables are not checked by C this routine. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WNLE C C ---------- C PARAMETERS C ---------- C INTEGER KEY, M REAL SCALR, X C DIMENSION X(M) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER I C C***FIRST EXECUTABLE STATEMENT WNLE C C ... QUICK CHECK (VECTORIZABLE) C KEY = 0 DO 10 I=1,M IF (X(I) .LE. SCALR) KEY = 1 10 CONTINUE C C ... IF CHECK FAILED THEN FIND INDEX OF FIRST FAILURE C IF (KEY .NE. 0) THEN KEY = 0 DO 20 I=1,M IF (X(I) .LE. SCALR) THEN KEY = I GO TO 30 ENDIF 20 CONTINUE ENDIF C 30 CONTINUE RETURN END SUBROUTINE WSCTR (N, Y, INDX, X) C***BEGIN PROLOGUE WSCTR C***PURPOSE Performs a scatter by index on a given vector C***LIBRARY VFNLIB C***TYPE SINGLE PRECISION (WSCTR-S, DWSCTR-D) C***AUTHOR SAUNDERS, B. V., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C BOISVERT, R. F., (NIST) C COMPUTING AND APPLIED MATHEMATICS LABORATORY C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C***DESCRIPTION C C WSCTR performs a scatter by index on a given vector based on the C indices provided in an array: X(INDX(i)) = Y(i), i=1,...,N. C C C P A R A M E T E R S C C N (Input) Integer (N .GT. 0) C The number of elements in the input vector Y. C C Y (Input) Real array of length N C The compressed vector. C C INDX (Input) Integer array of length N C Array of indices of positions elements of compressed C vector will occupy when scattered into vector X. C C X (Output) Real array C The vector to receive the scattered elements of Y. C X(INDX(i)) = Y(i), i=1,...,N. Only elements of X C whose indices are in INDX are changed. C C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 910226 DATE WRITTEN C***END PROLOGUE WSCTR C C C ---------- C PARAMETERS C ---------- C INTEGER INDX, N REAL Y, X C DIMENSION Y(N), X(*), INDX(N) C C --------------- C LOCAL VARIABLES C --------------- C INTEGER J C C***FIRST EXECUTABLE STATEMENT WGTHR C IF (N .LE. 0) RETURN C DO 10 J=1,N X(INDX(J)) = Y(J) 10 CONTINUE C RETURN END