C IQPACK - FORTRAN SUBROUTINES FOR THE WEIGHTS OF C INTERPOLATORY QUADRATURES C C FOR A DETAILED DESCRIPTION OF THESE ROUTINES SEE THE PAPER C WITH THE ABOVE TITLE - C C GIVEN A SET OF DISTINCT KNOTS, T, AND THEIR MULTIPLICITIES MLT, C THIS PACKAGE COMPUTES THE WEIGHTS D OF THE INTERPOLATORY C J,I C QUADRATURE FORMULA C C (I) C SUM SUM D F (T(J)), C J=1,NT I=0,MLT(J)-1 J,I C C (I) C WHERE F IS THE I-TH DERIVATIVE OF F, WHERE THE QUADRATURE C IS TO APPROXIMATE C C INTEGRAL F(T)W(T) DT, C [A,B] C C AND WHERE W(T) IS A WEIGHT FUNCTION. FOR CERTAIN CLASSICAL WEIGHT C FUNCTIONS, LISTED BELOW, NO OTHER INFORMATION IS NEEDED. HOWEVER C THE PACKAGE CAN COMPUTE THE QUADRATURE WEIGHTS CORRESPONDING TO C ANY W(T) FOR WHICH THE ZERO-TH MOMENT AND THE (TRIDIAGONAL C SYMMETRIC) JACOBI MATRIX ASSOCIATED WITH THE POLYNOMIALS C ORTHOGONAL ON [A,B] WITH RESPECT TO W(T), ARE KNOWN. (A UTILITY C ROUTINE IS SUPPLIED TO PROVIDE THIS INFORMATION FOR CLASSICAL C WEIGHT FUNCTIONS). KNOTS AND WEIGHTS OF GAUSS QUADRATURES C WITH NO MULTIPLE KNOTS CAN ALSO BE COMPUTED. C C THE PACKAGE IS AN IMPLEMENTATION OF THE METHOD DESCRIBED IN C C "CALCULATION OF THE WEIGHTS OF INTERPOLATORY QUADRATURES", C J. KAUTSKY AND S. ELHAY, NUMER MATH 40 (1982) 407-422, C C TOGETHER WITH VARIOUS UTILITY ROUTINES. WEIGHTS TO SOME OR ALL THE C KNOTS CAN BE COMPUTED. C C TABLE OF CLASSICAL WEIGHT FUNCTIONS C C KIND INTERVAL WEIGHT FUNCTION NAME C 1 (A,B) ONE LEGENDRE C 2 (A,B) ((B-X)*(X-A))**(-HALF) Chebyshev Type 1 C 3 (A,B) ((B-X)*(X-A))**ALPHA GEGENBAUER C 4 (A,B) (B-X)**ALPHA*(X-A)**BETA JACOBI C 5 (A,INF) (X-A)**ALPHA*EXP(-B*(X-A)) GEN LAGUERRE C 6 (-INF,INF) ABS(X-A)**ALFA*EXP(-B*(X-A)**2) GEN HERMITE C 7 (A,B) ABS(X-(A+B)/TWO)**ALFA EXPONENTIAL C 8 (A,INF) (X-A)**ALFA*(B+X)**BETA RATIONAL C C THE VALUES B=1 AND C A=-1 FOR WEIGHT FUNCTIONS 1,2,3,4,7 C A= 0 FOR WEIGHT FUNCTIONS 5,6,8 C WILL BE REFERRED TO AS THE DEFAULT VALUES. C C WE ALSO DEFINE DEL AS C (A+B)/2 FOR WEIGHT FUNCTIONS 1,2,3,4,7 C A FOR WEIGHT FUNCTIONS 5,6,8 C C IQPACK INDEX C ------------ C C LEGEND C ------ C GENERALLY I = THIS QUANTITY IS INPUT TO THIS ROUTINE C O = THIS QUANTITY IS OUTPUT FROM THIS ROUTINE C KNOTS - M = MULTIPLE KNOTS ALLOWED C S = ONLY SIMPLE KNOTS ALLOWED C WEIGHTS - C = COMPUTED C QF - I = ANY INTERPOLATORY QUADRATURE FORMULA C G = GAUSSIAN QUADRATURE FORMULA C EVAL - Y = THE QUADRATURE SUM IS FORMED C N = THE QUADRATURE SUM IS NOT FORMED C PRINT - Y = THE KNOTS AND WEIGHTS OF THE QUADRATURE FORMULA ARE C OPTIONALLY PRINTED AND A CHECK OF THE MOMENTS IS C OPTIONALLY PRINTED C N = NO PRINTING POSSIBLE C A,B - A = ANY VALID VALUES OF THE WEIGHT FUNCTION PARAMETERS A,B C ALLOWED C D = ONLY THE DEFAULT VALUES OF A,B ALLOWED C C USER ROUTINES C ------------- C NAME KNOTS WEIGHTS QF EVAL PRINT A,B C ------ C CEGQF SO C G Y N A C CEGQFS SO C G Y N D C CGQF SO OC G N Y A C CGQFS SO OC G N Y D C CDGQF SO OC G N N D C SGQF SO OC G N N - C CLIQF SI OC I N Y A C CLIQFS SI OC I N Y D C CEIQF MI C I Y N A C CEIQFS MI C I Y N D C CIQF MI CO I N Y A C CIQFS MI CO I N Y D C EIQF MI I I Y N - C EIQFS SI I I Y N - C CAWIQ MI C I N N D C C UTILITY AND AUXILLIARY ROUTINES C ------------------------------- C CLASS COMPUTE THE ZERO-TH MOMENT AND JACOBI MATRIX FOR C A CLASSICAL WEIGHT FUNCTION C WM COMPUTE THE MOMENTS OF A CLASSICAL WEIGHT FUNCTION C PARCHK CHECK THAT THE PARAMETER VALUES ARE VALID FOR THIS C WEIGHT FUNCTION C CHKQFS CHECK AND OPTIONALLY PRINT A MOMENTS CHECK OF A QF C AND OPTIONALLY PRINT THE KNOTS AND WEIGHTS. DEFAULT C VALUES OF A,B ONLY C CHKQF CHECK AND OPTIONALLY PRINT A MOMENTS CHECK OF A QF C AND OPTIONALLY PRINT THE KNOTS AND WEIGHTS. ANY C VALID VALUES OF A,B ALLOWED C SCT SCALE THE KNOTS OF A QF FOR ANY VALID A,B TO THOSE C FOR THE DEFAULT VALUES OF A,B C SCQF SCALE A CLASSICAL WEIGHT FUNCTION QF WITH DEFAULT C VALUES FOR A,B TO THOSE FOR ANY VALID A,B C SCMM SCALE THE MOMENTS OF A CLASSICAL WEIGHT FUNCTION C WITH DEFAULT VALUES FOR A,B TO THOSE FOR ANY VALID C A,B C WTFN COMPUTE THE VALUES OF A CLASSICAL WEIGHT FUNCTION C AT GIVEN POINTS C CWIQD FIND ALL THE WEIGHTS TO 1 MULTIPLE KNOT OF A QF C IMTQLX ORTHOGONALLY DIAGONALIZE A JACOBI MATRIX C MACHEP COMPUTE MACHINE EPSILON C DGAMMA COMPUTE DOUBLE PRECISION GAMMA FUNCTION C C THE FOLLOWING IS A LIST OF PARAMETERS USED THROUGHOUT THE PACKAGE C WHICH ALWAYS HAVE THE SAME MEANING. C C NT NUMBER OF DISTINCT KNOTS. MUST BE .GE.1. C T KNOT ARRAY. C MLT MULTIPLICITY ARRAY. T(J) HAS MULTIPLICITY MLT(J). C NWTS DIMENSION OF THE ARRAY CONTAINING THE WEIGHTS. C WTS ARRAY CONTAINING THE WEIGHTS. C NDX FLAGS AND POINTERS ARRAY. THE PACKAGE HAS BEEN DESIGNED TO C (1) TREAT ALL OR ONLY SOME OF THE KNOTS SUPPLIED AS INCLUDED IN C THE QUADRATURE, C (2) COMPUTE THE WEIGHTS FOR ALL OR ONLY SOME OF THE KNOTS C INCLUDED IN THE QUADRATURE, C (3) TO PACK THE WEIGHTS IN THE OUTPUT ARRAY IN VARIOUS (POSSIBLY C FOUR) DIFFERENT WAYS. C NDX INDICATES THE STATUS OF EACH KNOT AND POINTS TO THE LOCATION C OF THAT KNOT IN THE WTS ARRAY. ITS USE IS DESCRIBED IN CAWIQ. C IN MOST STRAIGHTFORWARD APPLICATIONS THE USER WILL ONLY NEED TO C DIMENSION THE ARRAY. THE PACKAGE WILL DO THE REST. C KEY WEIGHTS ARRAY STRUCTURE FLAG. WILL USUALLY BE SET 1. USE C DESCRIBED IN CAWIQ. C KIND AN INTEGER 0.LE.KIND.LE.8 SPECIFYING WHICH WEIGHT FUNCTION C IS TO BE USED. KIND=0 INDICATES THAT THE WEIGHT FUNCTION IS OF A C TYPE NOT LISTED IN THE TABLE BELOW OF CLASSICAL WEIGHT C FUNCTIONS. FOR KIND=0 THE USER MUST SUPPLY THE C JACOBI MATRIX AND ANY MOMENTS WHICH ARE REQUIRED. C ALPHA C BETA C A C B THE WEIGHT FUNCTION AND/OR INTERVAL PARAMETERS. ANY ONE MAY C BE REPLACED BY A DUMMY VARIABLE IF THE WEIGHT FUNCTION IS C INDEPENDENT OF IT. C NWF AN INTEGER SPECIFYING THE DIMENSION OF THE WORKFIELD WF. C MINIMUM VALUES FOR NWF ARE GIVEN IN THE DESCRIPTION OF EACH C ROUTINE THAT USES A WORKFIELD. C WF FLOATING POINT WORKFIELD ARRAY TO BE SUPPLIED BY THE USER. C NIWF AN INTEGER SPECIFYING THE DIMENSION OF IWF C IWF INTEGER TYPE WORKFIELD ARRAY TO BE SUPPLIED BY THE USER. C QFSUM VARIABLE RETURNING THE VALUE OF THE QUADRATURE SUM. C F A USER SUPPLIED FUNCTION INVOKED BY A STATEMENT LIKE Y=F(X,I). C IT RETURNS THE VALUE OF THE I-TH DERIVATIVE OF F AT X (ZERO-TH C DERIVATIVE=FUNCTION). THE FUNCTION SHOULD BE CAPABLE C OF RETURNING DERIVATIVES OF ALL ORDERS UP TO MMAX-1 WHERE C MMAX IS THE MAXIMUM MULTIPLICITY USED AT THE KNOTS. THE ACTUAL C PARAMETER USED IN THE CALL TO ROUTINE EIQF, EIQFS, CEIQF AND C CEIQFS MUST BE DECLARED IN AN EXTERNAL STATEMENT IN THE CALLING C PROGRAM C LO INTEGER VARIABLE USED TO CONTROL OUTPUT. IF LO IS SET TO ZERO C THEN THERE WILL BE NO OUTPUT PRINTED. IF LO IS NON-ZERO THEN C ABS(LO) WILL BE THE LOGICAL UNIT NUMBER TO WHICH ALL OUTPUT C IS DIRECTED. WHEN LO IS NEGATIVE WEIGHTS ONLY WILL BE PRINTED C AND WHEN LO IS POSITIVE THE WEIGHTS AND A CHECK OF THE MOMENTS C WILL BE PRINTED. IN SOME ROUTINES LO.EQ.0 WILL CAUSE A MOMENTS C CHECK TO BE COMPUTED EVEN THOUGH THERE IS NO PRINT WHILE IN C OTHERS LO.EQ.0 WILL CAUSE ONLY THE WEIGHTS TO BE COMPUTED. SEE C INDIVIDUAL ROUTINES FOR DETAILS. C C THROUGHOUT THE COMMENTS IN THIS PACKAGE C N...IS THE NUMBER OF KNOTS COUNTED ACCORDING TO THEIR C MULTIPLICITIES, C MMAX...MAXIMUM OF THE MLT(J) C RMAX...MAXIMUM OF 2*MMAX AND N+1 C NSTAR...INTEGER PART OF (N+1)/2 C C ERROR CONDITIONS ARE INDICATED BY THE VARIABLE IER BEING C RETURNED WITH A NON-ZERO VALUE. C C IER = 1 ALPHA.GT.-1 FALSE C 2 FOR KIND.LT.8 BETA.GT.-1 IS FALSE C 3 FOR KIND=8 NEED BETA.LT.(ALPHA+BETA+2*N).LT.0 C TO COMPUTE N ELEMENTS OF THE JACOBI MATRIX. C 4 UNKNOWN WEIGHT FUNCTION. CANNOT GENERATE C JACOBI MATRIX C 5 GAMMA FUNCTION AND MACHINE PARAMETERS ARE NOT C MATCHED IN ACCURACY C 6 ZERO LENGTH INTERVAL (KIND=1,2,3,4,7) C 7 B.LE.0 FOR KIND=5,6 C 8 A+B.LE.0 FOR KIND=8 C 9 NOT ENOUGH INTEGER WORKFIELD. NIWF=2*NT WILL DO C 10 DIMENSION OF WEIGHTS ARRAY TOO SMALL C 11 JACOBI MATRIX NOT DIAGONALIZED SUCCESSFULLY C 12 SIZE OF JACOBI MATRIX TOO SMALL FOR NUMBER OF C WEIGHTS C 13 ZERO-TH MOMENT OF WEIGHTS FUNCTION IS NOT > 0 C 14 KNOTS NOT DISTINCT C 15 SOME KNOT HAS MULTIPLICITY < 1 C 16 POINTERS FOR WGHTS ARRAY CONTRADICTORY C 17 0 < ABS(KEY) < 5 FALSE (SEE CAWIQ OR EIQF) C 18 NUMBER OF KNOTS < 1 C -K,K>0 AT LEAST K LOCATIONS ARE REQUIRED IN THE C FLOATING-POINT WORKFIELD IN ORDER TO COMPLETE C THE CURRENT TASK. C C SUBROUTINES AND THEIR CALL SEQUENCES C C CALL CEGQFS(NT,KIND,ALPHA,BETA,F,QFSUM, C 1 NWF,WF,NIWF,IWF,IER) C CALL CEGQF(NT,KIND,ALPHA,BETA,A,B,F,QFSUM, C 1 NWF,WF,NIWF,IWF,IER) C CALL CGQF(NT,T,WTS,KIND,ALPHA,BETA,A,B,LO, C 1 NWF,WF,NIWF,IWF,IER) C CALL CGQFS(NT,T,WTS,KIND,ALPHA,BETA,LO, C 1 NWF,WF,NIWF,IWF,IER) C CALL CDGQF(NT,T,WTS,KIND,ALPHA,BETA,NWF,WF,IER) C CALL SGQF(NT,T,WTS,AJ,BJ,ZEMU,IER) C CALL CLIQFS(NT,T,WTS,KIND,ALPHA,BETA, C 1 LO,NWF,WF,NIWF,IWF,IER) C CALL CLIQF(NT,T,WTS,KIND,ALPHA,BETA,A,B, C 1 LO,NWF,WF,NIWF,IWF,IER) C CALL CEIQFS(NT,T,MLT,KIND,ALPHA,BETA,F,QFSUM C 1 ,NWF,WF,NIWF,IWF,IER) C CALL CEIQF(NT,T,MLT,KIND,ALPHA,BETA,A,B,F,QFSUM C 1 ,NWF,WF,NIWF,IWF,IER) C CALL CIQFS(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA C 1 ,LO,NWF,WF,IER) C CALL CIQF(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA,A,B C 1 ,LO,NWF,WF,IER) C CALL EIQF(NT,T,MLT,WTS,NWTS,NDX,KEY,F,QFSUM,IER) C CALL EIQFS(NT,T,WTS,F,QFSUM,IER) C CALL CAWIQ(NT,T,MLT,NWTS,WTS,NDX,KEY C 1 ,NST,AJ,BJ,JDF,ZEMU,NWF,WF,IER) C CALL CWIQD(M,NM,L,V,XK,NSTAR,PHI,A,WF,Y,R,Z,D) C CALL CLASS(KIND,M,ALPHA,BETA,BJ,AJ,ZEMU,IER) C CALL WM(W,M,KIND,ALPHA,BETA,IER) C CALL PARCHK(KIND,M,ALPHA,BETA,IER) C CALL CHKQFS(T,WTS,MLT,NT,NWTS,NDX,KEY,W,MOP,MEX, C 1 KIND,ALPHA,BETA,LO,E,ER,QM,IER) C CALL CHKQF(T,WTS,MLT,NT,NWTS,NDX,KEY,WF,MOP,MEX,KIND, C 1 ALPHA,BETA,LO,E,ER,QM,NWF,A,B,IER) C CALL SCT(NT,T,ST,KIND,A,B,IER) C CALL SCQF(NT,T,MLT,WTS,NWTS,NDX,SWTS,ST, C 1 KIND,ALPHA,BETA,A,B,IER) C CALL SCMM(W,M,KIND,ALPHA,BETA,A,B,IER) C CALL WTFN(T,W,NT,KIND,ALPHA,BETA,IER) C CALL IMTQLX(N,D,E,Z,IER) C CALL MACHEP (X) C Y=DGAMMA(X) C C---------------------------------------------------------------------- C C IN THE DESCRIPTIONS OF THE ROUTINES BELOW ALL C THE INPUT AND OUTPUT PARAMETERS ARE INDICATED BY C THE SINGLE LETTER I OR O ALIGNED TO EACH VARIABLE IN THE C CALLING SEQUENCE. A * INDICATES THAT THE VARIABLE IS C SOMETIMES SET ON INPUT AND SOMETIMES SET BY THE ROUTINE. C SUBROUTINE CEGQF(NT,KIND,ALPHA,BETA,A,B,F,QFSUM 1,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 c C ROUTINE TO: C 1. COMPUTE ALL THE KNOTS AND WEIGHTS OF CLASSICAL WEIGHT C FUNCTION GAUSS QUADRATURE FORMULA WITH ALL SIMPLE KNOTS C FOR ANY VALID VALUES OF A AND B C 2. EVALUATE THE QUADRATURE SUM C C INPUT AND OUTPUT VARIABLES - C C I I I I I I I O C SUBROUTINE CEGQF(NT,KIND,ALPHA,BETA,A,B,F,QFSUM C 1,NWF,WF,NIWF,IWF,IER) C I O I O O C C THE USER SUPPLIES A FUNCTION F, WHICH MUST BE DECLARED IN AN C EXTERNAL STATEMENT IN THE CALLING PROGRAM, AND WHICH RETURNS C VALUES OF F. C C NEED NWF .GE. 2*NT C NIWF .GE. 2*NT C C FUNCTIONS AND SUBROUTINES REFERENCED - CGQF EIQFS F DOUBLE PRECISION A,ALPHA,B,BETA,QFSUM,WF,F INTEGER IER,KIND,LU,NA,NB,NC,NIWF,NT,NWF,IWF DIMENSION WF(NWF),IWF(NIWF) EXTERNAL F IER=0 IF(NIWF.LT.2*NT) THEN IER=9 RETURN ENDIF IF(NWF.LT.2*NT) THEN IER=-2*NT RETURN ENDIF C SET WORKFIELD FOR WEIGHTS AND KNOTS 10 LU=0 NA=1 NB=NA+NT NC=NB+NT+1 CALL CGQF(NT,WF(NB),WF(NA),KIND,ALPHA,BETA, 1 A,B,LU,NWF-NC,WF(NC),NIWF,IWF,IER) IF(IER.NE.0) RETURN C EVALUATE THE QUADRATURE SUM CALL EIQFS(NT,WF(NB),WF(NA),F,QFSUM,IER) RETURN END SUBROUTINE CEGQFS(NT,KIND,ALPHA,BETA,F,QFSUM 1,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 C C ROUTINE TO: C 1. COMPUTE ALL THE KNOTS AND WEIGHTS OF CLASSICAL WEIGHT C FUNCTION GAUSS QUADRATURE FORMULA WITH ALL SIMPLE KNOTS C FOR THE DEFAULT VALUES OF A AND B C 2. EVALUATE THE QUADRATURE SUM C C INPUT AND OUTPUT VARIABLES - C C I I I I I O C SUBROUTINE CEGQFS(NT,KIND,ALPHA,BETA,F,QFSUM C 1,NWF,WF,NIWF,IWF,IER) C I O I O O C C F MUST BE DECLARED IN AN EXTERNAL STATEMENT C IN THE CALLING PROGRAM. C C NEED NWF .GE. 2*NT C NIWF .GE. 2*NT C C FUNCTIONS AND SUBROUTINES REFERENCED - CGQFS EIQFS F DOUBLE PRECISION ALPHA,BETA,QFSUM,WF,F INTEGER IER,KIND,LU,NA,NB,NC,NIWF,NT,NWF,IWF DIMENSION WF(NWF),IWF(NIWF) EXTERNAL F C CHECK THERE IS ENOUGH FLOATING POINT AND INTEGER WORKSPACE IER=0 IF(NIWF.LT.2*NT) THEN IER=9 RETURN ENDIF IF(NWF.LT.2*NT) THEN IER=-2*NT RETURN ENDIF C ASSIGN WORKSPACE FOR KNOTS AND WEIGHTS 10 LU=0 NA=1 NB=NA+NT NC=NB+NT+1 CALL CGQFS(NT,WF(NB),WF(NA),KIND,ALPHA,BETA,LU, 1 NWF-NC,WF(NC),NIWF,IWF,IER) IF(IER.NE.0) RETURN C EVALUATE THE QUADRATURE SUM CALL EIQFS(NT,WF(NB),WF(NA),F,QFSUM,IER) RETURN END SUBROUTINE CGQF(NT,T,WTS,KIND,ALPHA,BETA,A,B,LO, 1NWF,WF,NIWF,IWF,IER) c*********************************************************************72 C C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF A GAUSS QF WITH C 1. A CLASSICAL WEIGHT FUNCTION WITH ANY VALID A,B C 2. ONLY SIMPLE KNOTS C 3. OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE AND PRINT KNOTS AND WEIGHTS. PRINT MOMENTS CHECK C LO.EQ.0...COMPUTE KNOTS AND WEIGHTS. PRINT NOTHING C LO.LT.0...COMPUTE AND PRINT KNOTS AND WEIGHTS. NO MOMENTS CHECK. C C INPUT AND OUTPUT VARIABLES - C I O O I I I I I I C SUBROUTINE CGQF(NT,T,WTS,KIND,ALPHA,BETA,A,B,LO, C I O I O O C 1NWF,WF,NIWF,IWF,IER) C C NEED NWF.GE. (9*NT+13) IF LO.GT.0 C (2*NT) IF LO.EQ.0 C (3*NT+4) IF LO.LT.0 C IWF...DIMENSION MUST BE .GE. 2*NT C C USE ROUTINE EIQFS TO EVALUATE A QUADRATURE COMPUTED BY CGQF. C C FUNCTIONS AND SUBROUTINES REFERENCED - CDGQF CHKQF SCQF DOUBLE PRECISION A,ALPHA,B,BETA,T,WF,WTS INTEGER I,IER,KEY,KIND,LEX,LO,M,MEX,MMEX,MOP,NAI,NBI,NE,NER INTEGER NILAST,NIWF,NLAST,NQM,NT,NW,NWF,IWF DIMENSION T(NT),WTS(NT),WF(NWF),IWF(NIWF) C C CHECK THERE IS ENOUGH WORKFIELD AND ASSIGN WORKFIELD IER=0 KEY=1 MOP=2*NT M=MOP+1 MEX=M+2 MMEX=MAX(MEX,1) LEX=MOP IF(LO.NE.0)LEX=MEX+NT+1 IF(LO.LE.0)MEX=0 NE=1 NER=NE+MEX NQM=NER+MEX NW=NQM+MEX NLAST=NW-1 LEX=LEX+3*MEX C EXIT IF INSUFFICIENT WORKFIELD IF(NIWF.LT.2*NT)IER=9 IF(NWF.LT.LEX)IER=-LEX IF(IER.NE.0) RETURN C C COMPUTE THE GAUSS QF FOR DEFAULT VALUES OF A,B CALL CDGQF(NT,T,WTS,KIND,ALPHA,BETA, 1NWF,WF,IER) C EXIT IF ERROR IF(IER.NE.0) RETURN C C PREPARE TO SCALE QF TO OTHER WEIGHT FUNCTION WITH VALID A,B C SET UP INTEGER WORK FIELDS NAI=1 NBI=NAI+NT NILAST=NBI+NT-1 DO 10 I=1,NT IWF(NAI+I-1)=1 IWF(NBI+I-1)=I 10 CONTINUE C IWF(NAI) IS THE MLT ARRAY. ALL KNOTS MULT=1 C IWF(NBI) IS THE NDX ARRAY. NDX(I)=I C SCALE THE QUADRATURE CALL SCQF(NT,T,IWF(NAI),WTS,NT,IWF(NBI),WTS,T, 1 KIND,ALPHA,BETA,A,B,IER) C C EXIT IF ERROR OR IF NO PRINT REQUIRED IF(IER.NE.0.OR.LO.EQ.0) RETURN C CALL CHKQF(T,WTS,IWF(NAI),NT,NT,IWF(NBI),KEY,WF(NW),MOP,MMEX, 1 KIND,ALPHA,BETA,LO,WF(NE),WF(NER),WF(NQM),NWF-NW,A,B,IER) RETURN END SUBROUTINE CGQFS(NT,T,WTS,KIND,ALPHA,BETA,LO, 1NWF,WF,NIWF,IWF,IER) c*********************************************************************72 C C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF A GAUSS QF WITH C 1. A CLASSICAL WEIGHT FUNCTION WITH DEFAULT VALUES FOR A,B C 2. ONLY SIMPLE KNOTS C 3. OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE AND PRINT KNOTS AND WEIGHTS. PRINT MOMENTS CHECK C LO.EQ.0...COMPUTE KNOTS AND WEIGHTS. PRINT NOTHING C LO.LT.0...COMPUTE AND PRINT KNOTS AND WEIGHTS. NO MOMENTS CHECK. C C INPUT AND OUTPUT VARIABLES - C I O O I I I I C SUBROUTINE CGQFS(NT,T,WTS,KIND,ALPHA,BETA,LO, C 1NWF,WF,NIWF,IWF,IER) C I O I O O C C NEED NWF.GE. (9*NT+13) IF LO.GT.0 C (2*NT) IF LO.EQ.0 C (3*NT+4) IF LO.LT.0 C IWF...DIMENSION MUST BE .GE. 2*NT C C USE ROUTINE EIQFS TO EVALUATE A QUADRATURE COMPUTED BY CGQFS. C C FUNCTIONS AND SUBROUTINES REFERENCED - CDGQF CHKQFS DOUBLE PRECISION ALPHA,BETA,T,WF,WTS INTEGER I,IER,KEY,KIND,LEX,LO,M,MEX,MMEX,MOP,NAI,NBI,NE,NER INTEGER NILAST,NIWF,NLAST,NQM,NT,NW,NWF,IWF DIMENSION T(NT),WTS(NT),WF(NWF),IWF(NIWF) C C CHECK THERE IS ENOUGH WORKFIELD AND ASSIGN WORKFIELD IER=0 KEY=1 MOP=2*NT M=MOP+1 MEX=M+2 MMEX=MAX(MEX,1) LEX=MOP IF(LO.NE.0)LEX=MEX+NT+1 IF(LO.LE.0)MEX=0 NE=1 NER=NE+MEX NQM=NER+MEX NW=NQM+MEX NLAST=NW-1 LEX=LEX+3*MEX C EXIT IF INSUFFICIENT WORKFIELD IF(NIWF.LT.2*NT)IER=9 IF(NWF.LT.LEX)IER=-LEX IF(IER.NE.0) RETURN C C COMPUTE THE GAUSS QF CALL CDGQF(NT,T,WTS,KIND,ALPHA,BETA, 1NWF,WF,IER) C EXIT IF ERROR OR IF NO PRINT REQUIRED IF(IER.NE.0.OR.LO.EQ.0) RETURN C C SET UP INTEGER WORK FIELDS NAI=1 NBI=NAI+NT NILAST=NBI+NT-1 DO 10 I=1,NT IWF(NAI+I-1)=1 IWF(NBI+I-1)=I 10 CONTINUE C IWF(NAI) IS THE MLT ARRAY. ALL KNOTS MULT=1 C IWF(NBI) IS THE NDX ARRAY. NDX(I)=I C CALL CHKQFS(T,WTS,IWF(NAI),NT,NT,IWF(NBI),KEY,WF(NW),MOP,MMEX, 1 KIND,ALPHA,BETA,LO,WF(NE),WF(NER),WF(NQM),IER) RETURN END SUBROUTINE CDGQF(NT,T,WTS,KIND,ALPHA,BETA, 1NWF,WF,IER) c*********************************************************************72 C C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF A GAUSS QF WITH C 1. A CLASSICAL WEIGHT FUNCTION WITH DEFAULT VALUES FOR A,B C 2. ONLY SIMPLE KNOTS C NO MOMENTS CHECK OR PRINTING DONE. C C INPUT AND OUTPUT VARIABLES - C I O O I I I C SUBROUTINE CDGQF(NT,T,WTS,KIND,ALPHA,BETA, C 1NWF,WF,IER) C I O O C C NWF... MUST BE .GE. 2*NT C C USE ROUTINE EIQFS TO EVALUATE A QUADRATURE COMPUTED BY CGQFS. C C FUNCTIONS AND SUBROUTINES REFERENCED - CLASS PARCHK SGQF DOUBLE PRECISION ALPHA,BETA,ZEMU,T,WF,WTS INTEGER IER,KIND,LEX,NA,NB,NLAST,NT,NWF DIMENSION T(NT),WTS(NT),WF(NWF) CALL PARCHK(KIND,2*NT,ALPHA,BETA,IER) C SET UP ARRAYS FOR DIAGONAL AND SUB-DIAGONAL OF JACOBI MATRIX NA=1 NB=NA+NT NLAST=NB+NT-1 LEX=2*NT IF(NWF.LT.LEX)IER=-LEX IF(IER.NE.0) RETURN C GET JACOBI MATRIX AND ZERO-TH MOMENT 10 CALL CLASS(KIND,NT,ALPHA,BETA,WF(NB),WF(NA),ZEMU,IER) IF(IER.NE.0) RETURN CALL SGQF(NT,T,WTS,WF(NA),WF(NB),ZEMU,IER) RETURN END SUBROUTINE SGQF(NT,T,WTS,AJ,BJ,ZEMU,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF A GAUSS QUADRATURE C FORMULA (WITH SIMPLE KNOTS) FROM THE JACOBI MATRIX AND THE ZERO-TH C MOMENT OF THE WEIGHT FUNCTION, USING THE GOLUB-WELSCH TECHNIQUE C C INPUT AND OUTPUT VARIABLES - C I O O I I I O C SUBROUTINE SGQF(NT,T,WTS,AJ,BJ,ZEMU,IER) C C INPUT PARAMETERS C AJ...DIAGONAL OF JACOBI MATRIX C BJ...SUB-DIAGONAL OF JACOBI MATRIX ( IN BJ(1)..BJ(NT-1) ) C ZEMU...ZERO-TH MOMENT OF WEIGHT FUNCTION C C OUTPUT PARAMETERS C AT OUTPUT T AND WTS CONTAIN THE KNOTS AND WEIGHTS C THE ARRAY BJ IS OVERWRITTEN DURING EXECUTION C C FUNCTIONS AND SUBROUTINES REFERENCED - IMTQLX MACHEP SQRT DOUBLE PRECISION ZEMU,AJ,BJ,PREC,T,WTS DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,NT DIMENSION T(NT),WTS(NT),AJ(NT),BJ(NT) C PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) COMMON /CTRLR/ PREC(10) IER=0 C C COMPUTE MACHINE EPSILON FOR IMTQLX CALL MACHEP(PREC(1)) C C EXIT IF ZERO-TH MOMENT NOT POSITIVE IF(ZEMU.LE.ZERO)IER=13 IF(IER.NE.0) RETURN C SET UP VECTORS FOR IMTQLX DO 10 I=1,NT T(I)=AJ(I) WTS(I)=ZERO 10 CONTINUE WTS(1)=SQRT(ZEMU) C C DIAGONALIZE JACOBI MATRIX CALL IMTQLX(NT,T,BJ,WTS,IER) C C CHECK FOR ERROR RETURN FROM IMTQLX IF(IER.EQ.0) GOTO 20 IER=11 RETURN C 20 DO 30 I=1,NT WTS(I)=WTS(I)**2 30 CONTINUE RETURN END SUBROUTINE CLIQFS(NT,T,WTS,KIND,ALPHA,BETA, 1LO,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 C C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF AN INTERPOLATORY C QF WITH C 1. A CLASSICAL WEIGHT FUNCTION WITH DEFAULT VALUES FOR A,B C 2. ONLY SIMPLE KNOTS C 3. OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. PRINT MOMENTS CHECK. C LO.EQ.0...COMPUTE WEIGHTS. PRINT NOTHING. C LO.LT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. C C INPUT AND OUTPUT VARIABLES - C I I O I I I C SUBROUTINE CLIQFS(NT,T,WTS,KIND,ALPHA,BETA, C 1LO,NWF,WF,NIWF,IWF,IER) C I I O I O O C C NEED NWF .GE. (5*N+9)/2 IF LO.LE.0 C (9*N+25)/2 IF LO.GT.0 C NIWF .GE. 2*NT C C USE ROUTINE EIQFS TO EVALUATE A QUADRATURE COMPUTED BY CLIQFS. C C FUNCTIONS AND SUBROUTINES REFERENCED - CIQFS DOUBLE PRECISION ALPHA,BETA,T,WF,WTS INTEGER I,IER,KEY,KIND,LO,NA,NB,NIWF,NT,NW,NWF,IWF DIMENSION T(NT),WTS(NT),WF(NWF),IWF(NIWF) C IER=0 IF(NIWF.GE.2*NT) GOTO 10 IER=9 RETURN 10 KEY=1 C SET UP WORKFIELD FOR MLT,NDX NA=1 NB=NA+NT NW=NB+NT DO 20 I=1,NT IWF(I)=1 20 CONTINUE CALL CIQFS(NT,T,IWF(NA),NT,WTS,IWF(NB),KEY,KIND,ALPHA,BETA 1,LO,NWF,WF,IER) RETURN END SUBROUTINE CLIQF(NT,T,WTS,KIND,ALPHA,BETA,A,B, 1LO,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE ALL THE KNOTS AND WEIGHTS OF AN INTERPOLATORY C QF WITH C 1. ONLY SIMPLE KNOTS AND C 2. A CLASSICAL WEIGHT FUNCTION WITH ANY VALID A,B C 3. OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. PRINT MOMENTS CHECK. C LO.EQ.0...COMPUTE WEIGHTS. PRINT NOTHING. C LO.LT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. C C INPUT AND OUTPUT VARIABLES - C I I O I I I I I C SUBROUTINE CLIQF(NT,T,WTS,KIND,ALPHA,BETA,A,B, C 1LO,NWF,WF,NIWF,IWF,IER) C I I O I O O C C NEED NWF .GE. (5*N+9)/2 IF LO.LE.0 C (9*N+25)/2 IF LO.GT.0 C NIWF .GE. 2*NT C C USE ROUTINE EIQFS TO EVALUATE A QUADRATURE COMPUTED BY CLIQF. C C FUNCTIONS AND SUBROUTINES REFERENCED - CIQF DOUBLE PRECISION A,ALPHA,B,BETA,T,WF,WTS INTEGER I,IER,KEY,KIND,LO,NA,NB,NIWF,NT,NW,NWF,IWF DIMENSION T(NT),WTS(NT),WF(NWF),IWF(NIWF) C IER=0 IF(NIWF.GE.2*NT) GOTO 10 IER=9 RETURN 10 KEY=1 C SET UP WORKFIELD FOR MLT,NDX NA=1 NB=NA+NT NW=NB+NT DO 20 I=1,NT IWF(I)=1 20 CONTINUE CALL CIQF(NT,T,IWF(NA),NT,WTS,IWF(NB),KEY,KIND,ALPHA,BETA,A,B 1,LO,NWF,WF,IER) RETURN END SUBROUTINE CEIQFS(NT,T,MLT,KIND,ALPHA,BETA,F,QFSUM 1,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 c C ROUTINE TO: C 1. COMPUTE AN INTERPOLATORY QF FOR CLASSICAL C WEIGHT FUNCTION WITH DEFAULT VALUES FOR A,B C 2. EVALUATE THE QUADRATURE SUM C C INPUT AND OUTPUT VARIABLES - C I I I I I I I O C SUBROUTINE CEIQFS(NT,T,MLT,KIND,ALPHA,BETA,F,QFSUM C 1,NWF,WF,NIWF,IWF,IER) C I O I O O C C NEED NWF .GE. NSTAR+RMAX+NT+3*(N+1) C NIWF .GE. NT C C FUNCTION F, MUST BE DECLARED IN AN EXTERNAL STATEMENT C IN THE CALLING PROGRAM. C C FUNCTIONS AND SUBROUTINES REFERENCED - CIQFS EIQF F DOUBLE PRECISION ALPHA,BETA,QFSUM,T,WF,F INTEGER IER,IWF,J,KEY,KIND,L,LEX,LU,M,MLT,MTM,N,NA,NIWF INTEGER NST,NT,NW,NWF DIMENSION T(NT),MLT(NT),WF(NWF),IWF(NIWF) EXTERNAL F C IER=0 IF(NIWF.GE.NT) GOTO 10 IER=9 RETURN 10 LU=0 N=0 MTM=MLT(1) DO 20 J=1,NT MTM=MAX(MTM,MLT(J)) 20 N=N+MLT(J) M=N+1 NST=M/2 L=MIN(2*MTM,M) LEX=NST+3*M+L+NT IF(NWF.GE.LEX) GOTO 30 IER=-LEX RETURN C INDECIES FOR WTS,NDX,WF (RESP) 30 NA=1 NW=NA+N KEY=1 CALL CIQFS(NT,T,MLT,N,WF(NA),IWF,KEY,KIND,ALPHA,BETA 1,LU,NWF-NW,WF(NW),IER) IF(IER.NE.0) RETURN CALL EIQF(NT,T,MLT,WF(NA),N,IWF,KEY,F,QFSUM,IER) RETURN END SUBROUTINE CEIQF(NT,T,MLT,KIND,ALPHA,BETA,A,B,F,QFSUM 1,NWF,WF,NIWF,IWF,IER) c*********************************************************************72 c C ROUTINE TO: C 1. COMPUTE AN INTERPOLATORY QF WITH CLASSICAL C WEIGHT FUNCTION WITH ANY VALID A,B C 2. EVALUATE THE QUADRATURE SUM C C INPUT AND OUTPUT VARIABLES - C I I I I I I I I I O C SUBROUTINE CEIQF(NT,T,MLT,KIND,ALPHA,BETA,A,B,F,QFSUM C 1,NWF,WF,NIWF,IWF,IER) C I O I O O C C NEED NWF .GE. NSTAR+RMAX+NT+3*(N+1) C NIWF .GE. NT C C FUNCTION F, MUST BE DECLARED IN AN EXTERNAL STATEMENT C IN THE CALLING PROGRAM. C C FUNCTIONS AND SUBROUTINES REFERENCED - CIQF EIQF F DOUBLE PRECISION A,ALPHA,B,BETA,QFSUM,T,WF,F INTEGER IER,IWF,J,KEY,KIND,L,LEX,LU,M,MLT,MTM,N,NA,NIWF INTEGER NST,NT,NW,NWF DIMENSION T(NT),MLT(NT),WF(NWF),IWF(NIWF) EXTERNAL F C IER=0 IF(NIWF.GE.NT) GOTO 10 IER=9 RETURN 10 LU=0 N=0 MTM=MLT(1) DO 20 J=1,NT MTM=MAX(MTM,MLT(J)) 20 N=N+MLT(J) M=N+1 NST=M/2 L=MIN(2*MTM,M) LEX=NST+3*M+L+NT IF(NWF.GE.LEX) GOTO 30 IER=-LEX RETURN C INDECIES FOR WTS,WF 30 NA=1 NW=NA+N KEY=1 CALL CIQF(NT,T,MLT,N,WF(NA),IWF,KEY,KIND,ALPHA,BETA,A,B 1,LU,NWF-NW,WF(NW),IER) IF(IER.NE.0) RETURN CALL EIQF(NT,T,MLT,WF(NA),N,IWF,KEY,F,QFSUM,IER) RETURN END SUBROUTINE CIQFS(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA 1,LO,NWF,WF,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE SOME OR ALL THE WEIGHTS OF A C QF FOR A CLASSICAL WEIGHT FUNCTION WITH DEFAULT VALUES OF A,B C AND A GIVEN SET OF KNOTS AND MULTIPLICITIES. THE WEIGHTS MAY BE C PACKED INTO THE OUTPUT ARRAY WTS ACCORDING TO A USER-DEFINED C PATTERN OR SEQUENTIALLY. THE ROUTINE WILL ALSO C OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. PRINT MOMENTS CHECK. C LO.EQ.0...COMPUTE WEIGHTS. PRINT NOTHING. NO MOMENTS CHECK. C LO.LT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. NO MOMENTS CHECK. C C INPUT AND OUTPUT VARIABLES - C I I I I O * I I I I C SUBROUTINE CIQFS(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA C 1,LO,NWF,WF,IER) C I I O O C C NEED NWF.GE.NSTAR+RMAX+2*(N+1) IF NO MOMENTS CHECK REQUIRED C NSTAR+RMAX+2*(2*N+5)IF A MOMENTS CHECK REQUIRED C KEY...AN INTEGER VARIABLE INDICATING THE STRUCTURE OF THE WTS C ARRAY. IT WILL NORMALLY BE SET TO 1. FOR MORE DETAILS SEE C THE COMMENTS IN CAWIQ. C NDX...AN INTEGER ARRAY OF DIMENSION NT USED TO INDEX THE OUTPUT C ARRAY WTS. IF KEY=1 NDX NEED NOT BE PRESET. FOR MORE C DETAILS SEE THE COMMENTS IN CAWIQ. C C FUNCTIONS AND SUBROUTINES REFERENCED - CAWIQ CHKQFS CLASS DOUBLE PRECISION ALPHA,BETA,T,WF,WTS,ZEMU INTEGER IER,IFL,J,JDF,K,KEY,KIND,L,LEX,LO,M,MEX,MLT,MMEX INTEGER MTM,N,NA,NB,ND,NDX,NE,NF,NST,NT,NW,NWF,NWTS DIMENSION T(NT),MLT(NT),WTS(NWTS),WF(NWF),NDX(NT) C IER=0 JDF=0 N=0 MTM=MLT(1) L=ABS(KEY) DO 20 J=1,NT IF(L.EQ.1) GOTO 10 IF(ABS(NDX(J)).EQ.0) GOTO 20 10 MTM=MAX(MTM,MLT(J)) N=N+MLT(J) 20 CONTINUE C N KNOTS WHEN COUNTED ACCORDING TO MULT IF(NWTS.GE.N) GOTO 30 IER=10 RETURN 30 M=N+1 MEX=2+M NST=M/2 IFL=1 IF(LO.LE.0)IFL=0 L=MIN(2*MTM,M) K=MAX(M,3*(MEX)*IFL) LEX=NST+M+L+K IF(NWF.GE.LEX) GOTO 40 IER=-LEX RETURN C C SET UP WORK FIELD INDECIES FOR CLASS AND CAWIQ 40 NA=1 NB=NA+NST NW=NB+NST C C GET JACOBI MATRIX CALL CLASS(KIND,NST,ALPHA,BETA,WF(NB),WF(NA),ZEMU,IER) IF(IER.NE.0) RETURN C C CALL WEIGHTS ROUTINE CALL CAWIQ(NT,T,MLT,N,WTS,NDX,KEY 1,NST,WF(NA),WF(NB),JDF,ZEMU,NWF-NW,WF(NW),IER) C RETURN IF NO PRINTING OR CHECKING REQUIRED 50 IF(IER.NE.0.OR.LO.EQ.0) RETURN MMEX=MEX*IFL ND=1 NE=ND+MMEX NF=NE+MMEX NW=NF+MMEX C CALL CHECKING ROUTINE CALL CHKQFS(T,WTS,MLT,NT,N,NDX,KEY,WF(NW),M-1,MEX,KIND, 1 ALPHA,BETA,LO,WF(ND),WF(NE),WF(NF),IER) RETURN END SUBROUTINE CIQF(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA,A,B 1,LO,NWF,WF,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE SOME OR ALL THE WEIGHTS OF A C QF FOR A CLASSICAL WEIGHT FUNCTION WITH ANY VALID A,B AND C A GIVEN SET OF KNOTS AND MULTIPLICITIES. THE WEIGHTS MAY BE C PACKED INTO THE OUTPUT ARRAY WTS ACCORDING TO A USER-DEFINED C PATTERN OR SEQUENTIALLY. THE ROUTINE WILL ALSO C OPTIONALLY PRINT KNOTS AND WEIGHTS AND A CHECK OF THE MOMENTS C C LO.GT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. PRINT MOMENTS CHECK. C LO.EQ.0...COMPUTE WEIGHTS. PRINT NOTHING. NO MOMENTS CHECK. C LO.LT.0...COMPUTE WEIGHTS. PRINT WEIGHTS. NO MOMENTS CHECK. C C INPUT AND OUTPUT VARIABLES - C I I I I O * I I I I I I C SUBROUTINE CIQF(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA,A,B C 1,LO,NWF,WF,IER) C I I O O C C NEED NWF.GE.NSTAR+RMAX+2*(N+1) IF NO MOMENTS CHECK REQUIRED C NSTAR+RMAX+5*N+NT+13 IF A MOMENTS CHECK IS REQUIRED C KEY...AN INTEGER VARIABLE INDICATING THE STRUCTURE OF THE WTS C ARRAY. IT WILL NORMALLY BE SET TO 1. THIS WILL CAUSE THE C WEIGHTS TO BE PACKED SEQUENTIALLY IN ARRAY WTS. C FOR MORE DETAILS SEE THE COMMENTS IN CAWIQ. C NDX...AN INTEGER ARRAY OF DIMENSION NT USED TO INDEX THE OUTPUT C ARRAY WTS. IF KEY=1 NDX NEED NOT BE PRESET. FOR MORE C DETAILS SEE THE COMMENTS IN CAWIQ. C C FUNCTIONS AND SUBROUTINES REFERENCED - CHKQF CIQFS SCQF SCT DOUBLE PRECISION A,ALPHA,B,BETA,T,WF,WTS INTEGER IER,IFL,J,K,KEY,KIND,L,LEX,LO,LU,M,MEX,MMEX,MLT,MTM INTEGER ND,NDX,NE,NF,NST,NT,NW,NWF,NWTS DIMENSION T(NT),MLT(NT),WTS(NWTS),WF(NWF),NDX(NT) C IER=0 M=1 MTM=1 L=ABS(KEY) DO 20 J=1,NT IF(L.EQ.1) GOTO 10 IF(ABS(NDX(J)).EQ.0) GOTO 20 10 MTM=MAX(MTM,MLT(J)) M=M+MLT(J) 20 CONTINUE IF(NWTS+1.GE.M) GOTO 30 IER=10 RETURN 30 NST=M/2 MEX=2+M IFL=1 IF(LO.LE.0)IFL=0 L=MIN(2*MTM,M) K=MAX(M,3*MEX*IFL) LEX=NST+L+K+M+(MEX+NT)*IFL IF(NWF.GE.LEX) GOTO 40 IER=-LEX RETURN 40 CONTINUE NST=1 NF=NST+NT C SCALE THE KNOTS TO DEFAULT A,B CALL SCT(NT,T,WF(NST),KIND,A,B,IER) IF(IER.NE.0) RETURN LU=0 CALL CIQFS(NT,WF(NST),MLT,NWTS,WTS,NDX,KEY, 1KIND,ALPHA,BETA,LU,NWF-NF+1,WF(NF),IER) IF(IER.NE.0) RETURN C DON'T SCALE USER'S KNOTS - ONLY SCALE WEIGHTS CALL SCQF(NT,WF(NST),MLT,WTS,NWTS,NDX,WTS,WF(NST), 1KIND,ALPHA,BETA,A,B,IER) IF(IER.NE.0.OR.LO.EQ.0) RETURN MMEX=MEX*IFL ND=1 NE=ND+MMEX NF=NE+MMEX NW=NF+MMEX CALL CHKQF(T,WTS,MLT,NT,NWTS,NDX,KEY,WF(NW),M-1,MEX,KIND, 1 ALPHA,BETA,LO,WF(ND),WF(NE),WF(NF),NWF-NW,A,B,IER) RETURN END SUBROUTINE EIQF(NT,T,MLT,WTS,NWTS,NDX,KEY,F,QFSUM,IER) c*********************************************************************72 c C ROUTINE TO EVALUATE AN INTERPOLATORY QF WITH KNOTS, WEIGHTS C AND INTEGRAND SUPPLIED. C ALL KNOTS FOR WHICH NDX(I).NE.0 ARE USED. C INPUT AND OUTPUT VARIABLES - C I I I I I I I I O O C SUBROUTINE EIQF(NT,T,MLT,WTS,NWTS,NDX,KEY,F,QFSUM,IER) C C ************************************************************** C * C * F.......A FUNCTION WITH 2 PARAMETERS F(X,I) C * TO BE SUPPLIED BY THE USER. IT MUST RETURN THE I-TH C * DERIVATIVE OF F, THE FUNCTION BEING INTEGRATED, AT X. C * I MUST RANGE FROM 0 (FOR THE FUNCTION VALUES) UP TO C * (THE MAXIMUM VALUE IN MLT)-1. THIS FUNCTION WILL ONLY C * BE CALLED WITH F AND ITS DERIVATIVES AT THE KNOTS SUPPLIED C * SO IT CAN GENERATE THE VALUES FOR F BY TABLE LOOKUP. C * THIS CAN BE ACHIEVED BY REPLACING THE LINE C * QFSUM=QFSUM+WTS(L+I-1)*F(T(J),I-1)/P C * WITH C * QFSUM=QFSUM+WTS(L+I-1)*F(T,J,I-1)/P C * WHERE THE FUNCTION F HAS THE KNOTS ARRAY T AS A PARAMETER C * AND THE REQUIRED KNOT IS INDICATED BY THE INDEX J. F IS C * CALLED ONLY FROM THIS ROUTINE AND EIQFS. C * C ************************************************************** C FUNCTIONS AND SUBROUTINES REFERENCED - F DOUBLE PRECISION P,QFSUM,T,WTS,F INTEGER I,IER,J,KEY,L,MLT,NDX,NT,NWTS DOUBLE PRECISION ZERO,HALF,ONE,TWO DIMENSION T(NT),MLT(NT),WTS(NWTS),NDX(NT) EXTERNAL F PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C IER=0 L=ABS(KEY) IF(L.GE.1.AND.L.LE.4) GOTO 10 IER=17 RETURN 10 QFSUM=ZERO DO 30 J=1,NT L=ABS(NDX(J)) IF(L.EQ.0) GOTO 30 P=ONE DO 20 I=1,MLT(J) QFSUM=QFSUM+WTS(L+I-1)*F(T(J),I-1)/P IF(KEY.GT.0) GOTO 20 P=P*I 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE EIQFS(NT,T,WTS,F,QFSUM,IER) c*********************************************************************72 c C ROUTINE TO EVALUATE AN INTERPOLATORY QF WITH ALL KNOTS SIMPLE AND C ALL KNOTS INCLUDED IN THE QUADRATURE. THIS ROUTINE WILL BE USED C TYPICALLY AFTER CLIQF OR CLIQFS HAVE BEEN CALLED. C INPUT AND OUTPUT VARIABLES - C I I I I O O C SUBROUTINE EIQFS(NT,T,WTS,F,QFSUM,IER) C C ************************************************************** C * C * F.......A FUNCTION WITH 2 PARAMETERS F(X,I) TO BE SUPPLIED C * BY THE USER. F MUST RETURN THE VALUE OF F, C * THE FUNCTION BEING INTEGRATED, AT X. C * I MAY BE A DUMMY VARIABLE BUT IS INCLUDED TO MAKE THIS C * DEFINITION CONFORM WITH THAT FOR F ELSEWHERE. THIS FUNCTION C * WILL ONLY BE CALLED WITH F AND ITS DERIVATIVES AT THE KNOTS C * SUPPLIED SO IT CAN GENERATE THE VALUES FOR F BY TABLE LOOKUP. C * THIS CAN BE ACHIEVED BY REPLACING THE LINE C * QFSUM=QFSUM+WTS(J)*F(T(J),0) C * WITH C * QFSUM=QFSUM+WTS(J)*F(T,J,0) C * WHERE THE FUNCTION F HAS THE KNOTS ARRAY T AS A PARAMETER C * AND THE REQUIRED KNOT IS INDICATED BY THE INDEX J. F IS C * CALLED ONLY FROM THIS ROUTINE AND EIQF. C * C ************************************************************** C FUNCTIONS AND SUBROUTINES REFERENCED - F DOUBLE PRECISION QFSUM,T,WTS,F INTEGER IER,J,NT DOUBLE PRECISION ZERO,HALF,ONE,TWO DIMENSION T(NT),WTS(NT) EXTERNAL F PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C IER=0 QFSUM=ZERO DO 10 J=1,NT QFSUM=QFSUM+WTS(J)*F(T(J),0) 10 CONTINUE RETURN END SUBROUTINE CAWIQ(NT,T,MLT,NWTS,WTS,NDX,KEY 1,NST,AJ,BJ,JDF,ZEMU,NWF,WF,IER) c*********************************************************************72 C C THIS ROUTINE, GIVEN A SET OF DISTINCT KNOTS, T, THEIR C MULTIPLICITIES MLT, THE JACOBI MATRIX ASSOCIATED WITH THE C POLYNOMIALS ORTHOGONAL WITH RESPECT TO THE WEIGHT FUNCTION W(X), C AND THE ZERO-TH MOMENT OF W(X) COMPUTES THE WEIGHTS OF THE C QUADRATURE C C (I) C SUM SUM D F (T(J)) C J=1,NT I=0,MLT(J)-1 J,I C C WHICH IS TO APPROXIMATE C C INTEGRAL F(T)W(T) DT C [A,B] C C THE ROUTINE MAKES VARIOUS CHECKS, AS INDICATED BELOW, SETS UP C VARIOUS VECTORS AND, IF NECESSARY, CALLS FOR THE DIAGONALIZATION C OF THE JACOBI MATRIX THAT IS ASSOCIATED WITH THE POLYNOMIALS C ORTHOGONAL WITH RESPECT TO W(X) ON [A,B]. THEN FOR EACH KNOT, THE C WEIGHTS OF WHICH ARE REQUIRED, IT CALLS THE ROUTINE CWIQD WHICH C COMPUTES ALL THE WEIGHTS FOR ANY GIVEN KNOT. C C INPUT AND OUTPUT VARIABLES - C I I I I O * I C SUBROUTINE CAWIQ(NT,T,MLT,NWTS,WTS,NDX,KEY C 1,NST,AJ,BJ,JDF,ZEMU,NWF,WF,IER) C I I I I I I O O C C PARAMETERS MARKED WITH A * MAY BE CHANGED BY THE SUBROUTINE C C *NDX THIS ARRAY ASSOCIATES WITH EACH DISTINCT KNOT T(J), C AN INTEGER NDX(J) WHICH IS SUCH THAT THE WEIGHT TO THE C I-TH DERIV VALUE OF F AT THE J-TH KNOT, IS STORED IN C C WTS(ABS(NDX(J))+I) J=1,2,...,NT, C I=0,1,2,...,MLT(J)-1 C ALSO: C NDX > 0 MEANS WEIGHTS ARE WANTED FOR THIS KNOT C < 0 MEANS WEIGHTS NOT WANTED FOR THIS KNOT BUT THE C KNOT IS TO BE INCLUDED IN THE QUADRATURE C = 0 MEANS IGNORE THIS KNOT COMPLETELY C KEY SWITCH INDICATING STRUCTURE OF OUTPUT ARRAYS WTS C AND NDX. C C ABS(KEY)= 1 SET UP POINTERS IN NDX FOR ALL KNOTS C IN T ARRAY (ROUTINE CAWIQ DOES THIS). C THE CONTENTS OF NDX ARE NOT TESTED C ON INPUT AND WEIGHTS ARE PACKED C SEQUENTIALLY IN WTS AS INDICATED C ABOVE C 2 SET UP POINTERS ONLY FOR KNOTS WHICH C HAVE NDX.NE.0 ON INPUT. ALL KNOTS C WHICH HAVE A NON-ZERO FLAG ARE C ALLOCATED SPACE IN WTS C 3 SET UP POINTERS ONLY FOR KNOTS WHICH C HAVE NDX.GT.0 ON INPUT. SPACE IN WTS C ALLOCATED ONLY FOR KNOTS WITH C NDX > 0 C 4 NDX ASSUMED TO BE PRESET AS POINTER C ARRAY ON INPUT C C AND KEY>0 FOR WEIGHTS WTS(J) REQUIRED IN STANDARD FORM C KEY<0 FOR J!WTS(J) REQUIRED C C NST DIMENSION OF JACOBI MATRIX. NST SHOULD BE BETWEEN (N+1)/2 C AND N. THE USUAL CHOICE WILL BE (N+1)/2 C *JDF FLAG TO INDICATE WHETHER JACOBI MATRIX NEEDS C DIAGONALIZING OR NOT C JDF= 0 DIAGONALIZATION REQUIRED C 1 DIAGONALIZATION NOT REQUIRED C C *AJ,BJ IT IS ASSUMED ON INPUT THAT C IF JDF = 0 THEN AJ CONTAINS THE DIAGONAL OF THE C JACOBI MATRIX AND C BJ CONTAINS THE SUBDIAGONAL. C C NOTE THAT BJ(NST) IS ALSO C DEFINED BUT NOT USED. C C IF JDF = 1 AJ CONTAINS THE EIGENVALUES OF C THE JACOBI MATRIX AND C BJ CONTAINS THE SQUARES OF THE C ELEMENTS OF THE 1ST ROW OF U THE C ORTHOGONAL MATRIX DIAGONALIZING THE C T C JACOBI MATRIX AS U D U . C ZEMU ZERO-TH MOMENT OF THE WEIGHT FUNCTION W(X) C NWF DIMENSION OF WORK FIELD. MUST HAVE NWF.GE.NST+RMAX+N C *IER ERROR FLAG: CODED AS FOLLOWS C 10 NWTS TOO SMALL C 11 JACOBI MATRIX NOT DIAGONALIZED SUCCESSFULLY C 12 NST TOO SMALL FOR N C 13 ZEMU > 0 FALSE C 14 KNOTS NOT DISTINCT C 15 MLT(J)<1 FOR SOME J C 16 POINTERS FOR WTS CONTRADICTORY C 17 0 < ABS(KEY) < 5 FALSE C 18 NT < 1 C C FUNCTIONS AND SUBROUTINES REFERENCED - CWIQD IMTQLX MACHEP SIGN DOUBLE PRECISION AJ,BJ,P,PREC,T,TMP,WF,WTS,ZEMU INTEGER I,IER,IERRX,IP,IPM,J,JDF,JJ,JP,K,KEY,L,M,MLT,MNM,MTJ,MTM INTEGER N,NDX,NK,NM,NMAX,NR,NST,NT,NW,NWF,NWTS,NY,NZ DOUBLE PRECISION ZERO,HALF,ONE,TWO DIMENSION T(NT),MLT(NT),WTS(NWTS) 1,NDX(NT),AJ(NST),BJ(NST),WF(NWF) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) COMMON /CTRLR/ PREC(10) C C COMPUTE MACHINE EPSILON CALL MACHEP(PREC(1)) C C EXIT IF NT < 1 IER=18 IF(NT.LT.1) RETURN IER=0 C C CHECK FOR INDISTINCT KNOTS IF(NT.EQ.1) GOTO 20 K=NT-1 DO 10 I=1,K TMP=T(I) L=I+1 DO 10 J=L,NT IF( ABS(TMP-T(J)).GT.PREC(1)) GOTO 10 IER=14 RETURN 10 CONTINUE C C CHECK MULTIPLICITIES, C SET UP VARIOUS USEFUL PARAMETERS AND C SET UP OR CHECK POINTERS TO WTS ARRAY 20 L=ABS(KEY) IF(L.GE.1.AND.L.LE.4) GOTO 30 IER=17 RETURN 30 K=1 GOTO(40,60,60,100),L 40 DO 50 I=1,NT NDX(I)=K IF(MLT(I).LT.1) GOTO 90 K=K+MLT(I) 50 CONTINUE N=K-1 GOTO 140 60 N=0 DO 80 I=1,NT IF(NDX(I).EQ.0) GOTO 80 IF(MLT(I).LT.1) GOTO 90 N=N+MLT(I) IF(NDX(I).LT.0.AND.L.EQ.3) GOTO 80 70 NDX(I)= SIGN(K,NDX(I)) K=K+MLT(I) 80 CONTINUE IF(K.LE.NWTS+1) GOTO 140 IER=10 RETURN 90 CONTINUE IER=15 RETURN 100 CONTINUE DO 120 I=1,NT IP=ABS(NDX(I)) IF(IP.EQ.0) GOTO 120 IPM=IP+MLT(I) IF(IPM.GT.NWTS) GOTO 130 IF(I.EQ.NT) GOTO 140 L=I+1 DO 110 J=L,NT JP=ABS(NDX(J)) IF(JP.EQ.0) GOTO 110 IF(JP.LE.IPM.AND.IP.LE.JP+MLT(J)) GOTO 130 110 CONTINUE 120 CONTINUE GOTO 140 130 IER=16 RETURN C C GET MAXIMUM MULTIPLICITY TO SEE IF ENOUGH STORE C IS AVAILABLE 140 MTM=1 DO 150 I=1,NT IF(NDX(I).GT.0)MTM=MAX(MTM,MLT(I)) 150 CONTINUE C C SET UP WORK FIELDS AND TEST SOME PARAMETERS IF(NST.LT.(N+1)/2)IER=12 NMAX=N+NST+MIN(2*MTM,N+1) IF(ZEMU.LE.ZERO)IER=13 IF(NMAX.GT.NWF)IER=-NMAX IF(IER.NE.0) RETURN C C TREAT A QF WITH 1 SIMPLE KNOT FIRST. IF(N.GT.1) GOTO 180 DO 160 I=1,NT IF(NDX(I).GT.0) GOTO 170 160 CONTINUE C WEIGHT NOT WANTED! DO NOTHING. RETURN 170 WTS(ABS(NDX(I)))=ZEMU RETURN C SKIP DIAGONALIZATION IF ALREADY DONE 180 IF(JDF.NE.0) GOTO 230 NW=1 C C SET UNIT VECTOR IN WORK FIELD TO GET BACK 1ST ROW OF Q DO 190 I=1,NST K=NW+I-1 WF(K)=ZERO 190 CONTINUE WF(NW)=ONE 200 IERRX=0 C C DIAGONALIZE JACOBI MATRIX CALL IMTQLX(NST,AJ,BJ,WF(NW),IERRX) C C CHECK FOR ERROR IF(IERRX.EQ.0) GOTO 210 IER=11 RETURN C C SIGNAL JACOBI MATRIX NOW DIAGONALIZED SUCCESSFULLY AND SAVE C SQUARES OF 1ST ROW OF U IN SUBDIAGONAL ARRAY C 210 JDF=1 DO 220 I=1,NST L=I+NW-1 BJ(I)=WF(L)**2 220 CONTINUE C C FIND ALL THE WEIGHTS FOR EACH KNOT FLAGGED 230 CONTINUE DO 280 I=1,NT IF(NDX(I).LE.0) GOTO 280 M=MLT(I) NM=N-M MNM=MAX(NM,1) L=MIN(M,NM+1) C C SET UP INDECIES FOR WORK FIELDS NK=NW+NST NY=NK+MNM NR=NY+M NZ=NR+L C C SET UP K-HAT MATRIX (FOR CWIQD) WITH KNOTS ACCORDING TO C THEIR MULTIPLICITIES K=NK DO 250 J=1,NT IF(NDX(J).EQ.0) GOTO 250 IF(J.EQ.I) GOTO 250 MTJ=MLT(J) DO 240 JJ=1,MTJ WF(K)=T(J) K=K+1 240 CONTINUE 250 CONTINUE C C SET UP RIGHT PRINCIPAL VECTOR ARRAY FOR WEIGHTS ROUTINE WF(NR)=ONE/ZEMU DO 260 J=2,L K=NR+J-1 WF(K)=ZERO 260 CONTINUE C C PICK UP POINTER FOR THE LOCATION OF THE WEIGHTS TO C BE OUTPUT K=NDX(I) C C FIND ALL THE WEIGHTS FOR THIS KNOT C CALL CWIQD(M,MNM,L,T(I),WF(NK),NST,AJ,BJ, 1WF(NW),WF(NY),WF(NR),WF(NZ),WTS(K)) IF(KEY.LT.0) GOTO 280 C C DIVIDE BY FACTORIALS FOR WEIGHTS IN STANDARD FORM IF(M.LT.2) GOTO 280 TMP=ONE IP=M-1 DO 270 J=2,IP P=J TMP=TMP*P L=K+J WTS(L)=WTS(L)/TMP 270 CONTINUE 280 CONTINUE RETURN END SUBROUTINE CWIQD(M,NM,L,V,XK,NSTAR,PHI,A,WF,Y,R,Z,D) c*********************************************************************72 c C ROUTINE TO COMPUTE ALL THE WEIGHTS TO A GIVEN KNOT. C FOR DETAILS SEE: C KAUTSKY AND ELHAY "CALCULATION OF THE WEIGHTS OF INTERPOLATORY C QUADRATURES", NUMER MATH 40 (1982) 407-422. C C VARIABLES NAMES USED CORRESPOND CLOSELY WITH THOSE IN THE ABOVE C MENTIONED PAPER C INPUT AND OUTPUT VARIABLES - C I I I I I I I I O O O O O C SUBROUTINE CWIQD(M,NM,L,V,XK,NSTAR,PHI,A,WF,Y,R,Z,D) C C M MULTIPLICITY OF THE KNOT IN QUESTION C NM MAX(N-M,1). N=NUMBER OF KNOTS USED COUNTED C ACCORDING TO MULTIPLICITY C L MIN(M,N-M+1) C V THE KNOT IN QUESTION C XK ALL BUT THE LAST M ENTRIES IN THE DIAGONAL OF K-HAT C NSTAR DIMENSION OF THE JACOBI MATRIX C PHI THE EIGENVALUES OF THE JACOBI MATRIX J C A THE SQUARE OF THE 1ST ROW OF THE ORTHOGONAL C MATRIX DIAGONALIZING J C WF WORK FIELD C Y Y-HAT C R VECTOR USED TO COMPUTE THE RIGHT PRINCIPAL VECTORS C Z VECTOR USED TO COMPUTE THE LEFT PRINCIPAL VECTORS C D OUTPUT ARRAY FOR THE WEIGHTS C OTHER VARIABLES ARE FOR TEMPORARY USE ONLY C DOUBLE PRECISION SUM,TMP,V,A,D,PHI,R,WF,XK,Y,Z DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,J,JR,K,L,LAST,M,MINIL,NM,NSTAR DIMENSION XK(NM),PHI(NSTAR),A(NSTAR),WF(NSTAR),Y(M), 1R(L),Z(M),D(M) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C C COMPUTE PRODUCTS REQUIRED FOR Y-HAT DO 20 J=1,NSTAR WF(J)=A(J) IF(NM.LT.1) GOTO 20 DO 10 I=1,NM WF(J)=WF(J)*(PHI(J)-XK(I)) 10 CONTINUE 20 CONTINUE C COMPUTE Y-HAT DO 40 I=1,M SUM=ZERO DO 30 J=1,NSTAR SUM=SUM+WF(J) WF(J)=WF(J)*(PHI(J)-V) 30 CONTINUE Y(I)=SUM 40 CONTINUE C IF N=1 THE RIGHT PRINCIPAL VECTOR IS ALREADY IN R. IF(NM.EQ.0) GOTO 80 C OTHERWISE COMPUTE THE R-PRINCIPAL VECTOR OF GRADE M-1 DO 70 I=1,NM TMP=V-XK(I) IF(L.EQ.1) GOTO 60 LAST=MIN(L,I+1) DO 50 JR=2,LAST J=LAST-JR+2 R(J)=TMP*R(J)+R(J-1) 50 CONTINUE 60 R(1)=TMP*R(1) 70 CONTINUE C COMPUTE LEFT PRINCIPAL VECTOR(S) AND WEIGHT FOR HIGHEST DERIV C THE FOLLOWING STATEMENT CONTAINS THE ONLY DIVISION IN THIS C ROUTINE. ANY TEST FOR OVERFLOW SHOULD BE MADE AFTER IT. 80 Z(1)=ONE/R(1) D(M)=Y(M)*Z(1) IF(M.EQ.1) RETURN C COMPUTE L-PRINCIPAL VECTOR DO110I=2,M SUM=ZERO IF(L.EQ.1) GOTO 100 MINIL=MIN(I,L) DO90J=2,MINIL K=I-J+1 SUM=SUM+R(J)*Z(K) 90 CONTINUE 100 Z(I)=-SUM*Z(1) 110 CONTINUE C ACCUMULATE WEIGHTS DO130I=2,M SUM=ZERO DO120J=1,I K=M-I+J SUM=SUM+Z(J)*Y(K) 120 CONTINUE K=M-I+1 130 D(K)=SUM RETURN END SUBROUTINE CLASS(KIND,M,ALPHA,BETA,BJ,AJ,ZEMU,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE THE DIAGONAL (AJ) AND SUB-DIAGONAL (BJ) C ELEMENTS OF THE ORDER M (TRIDIAGONAL SYMMETRIC) JACOBI MATRIX C ASSOCIATED WITH THE POLYNOMIALS ORTHOGONAL WITH RESPECT TO C THE WEIGHT FUNCTION SPECIFIED BY KIND. C FOR WEIGHT FUNCTIONS 1-7 M ELEMENTS ARE DEFINED IN BJ EVEN C THOUGH ONLY M-1 ARE NEEDED. FOR WEIGHT FUNCTION 8 BJ(M) IS C SET TO ZERO. C THE ZERO-TH MOMENT OF THE WEIGHT FUNCTION IS RETURNED IN C ZEMU. C INPUT AND OUTPUT VARIABLES - C I I I I O O O O C SUBROUTINE CLASS(KIND,M,ALPHA,BETA,BJ,AJ,ZEMU,IER) C C ERROR CODES ARE: C IER=1,2,3 PARAMETER RANGES ARE WRONG C IER=4 WEIGHT FUNCTION OF UNKNOWN TYPE. CANNOT GENERATE JACOBI C MATRIX C IER=5 GAMMA FUNCTION DOES NOT MATCH MACHINE PARAMETERS C IN ACCURACY C C FUNCTIONS AND SUBROUTINES REFERENCED - DGAMMA MACHEP SQRT PARCHK DOUBLE PRECISION AJ,BJ,A2B2,AB,ABA,ABI,ABJ,ABTI,ALPHA DOUBLE PRECISION APONE,BETA,TEMP,ZEMU,DGAMMA DOUBLE PRECISION ZERO,HALF,ONE,TWO DOUBLE PRECISION THREE,FOUR DOUBLE PRECISION PI INTEGER I,IER,KIND,M DIMENSION AJ(M),BJ(M) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) PARAMETER (THREE=3.0D0,FOUR=4.0D0) PARAMETER (PI=3.14159 26535 89793 23846 26433 83279 50 D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) CS PARAMETER (THREE=3.0E0,FOUR=4.0E0) CS PARAMETER (PI=3.14159 26535 89793 23846 26433 83279 50 E0) C CALL MACHEP(TEMP) CALL PARCHK(KIND,2*M-1,ALPHA,BETA,IER) IF(ABS(DGAMMA(HALF)**2-PI).GT.5.0D2*TEMP)IER=5 IF(IER.NE.0) RETURN C C LEG,CHEB,GEG,JAC,LAG,HERM,EXP,RAT GO TO( 20, 50, 70, 90,110, 130, 10,150),KIND 10 AB=ALPHA GOTO 30 20 AB=ZERO 30 ZEMU=TWO/(AB+ONE) DO 40 I=1,M AJ(I)=ZERO ABI=I+AB*MOD(I,2) ABJ=2*I+AB 40 BJ(I)=ABI*ABI/(ABJ*ABJ-ONE) GOTO 170 50 ZEMU=PI DO 60 I=1,M AJ(I)=ZERO 60 BJ(I)=HALF BJ(1)= SQRT(HALF) RETURN 70 AB=ALPHA*TWO ZEMU=TWO**(AB+ONE)*DGAMMA(ALPHA+ONE)**2/DGAMMA(AB+TWO) AJ(1)=ZERO BJ(1)=ONE/(TWO*ALPHA+THREE) DO 80 I=2,M AJ(I)=ZERO 80 BJ(I)=I*(I+AB)/(FOUR*(I+ALPHA)**2-ONE) GOTO 170 90 AB=ALPHA+BETA ABI=TWO+AB ZEMU=TWO**(AB+ONE)*DGAMMA(ALPHA+ONE)*DGAMMA(BETA+ONE)/ 1 DGAMMA(ABI) AJ(1)=(BETA-ALPHA)/ABI BJ(1)=FOUR*(ONE+ALPHA)*(ONE+BETA)/((ABI+ONE)*ABI*ABI) A2B2=BETA*BETA-ALPHA*ALPHA DO 100 I=2,M ABI=TWO*I+AB AJ(I)=A2B2/((ABI-TWO)*ABI) ABI=ABI**2 100 BJ(I)=FOUR*I*(I+ALPHA)*(I+BETA)*(I+AB)/((ABI-ONE)*ABI) GOTO 170 110 ZEMU=DGAMMA(ALPHA+ONE) DO 120 I=1,M AJ(I)=TWO*I-ONE+ALPHA 120 BJ(I)=I*(I+ALPHA) GOTO 170 130 ZEMU=DGAMMA((ALPHA+ONE)/TWO) DO 140 I=1,M AJ(I)=ZERO 140 BJ(I)=(I+ALPHA*MOD(I,2))/TWO GOTO 170 150 AB=ALPHA+BETA ZEMU=DGAMMA(ALPHA+ONE)*DGAMMA(-(AB+ONE))/DGAMMA(-BETA) APONE=ALPHA+ONE ABA=AB*APONE AJ(1)=-APONE/(AB+TWO) BJ(1)=-AJ(1)*(BETA+ONE)/(AB+TWO)/(AB+THREE) DO 160 I=2,M ABTI=AB+TWO*I AJ(I)=ABA+TWO*(AB+I)*(I-1) 160 AJ(I)=-AJ(I)/ABTI/(ABTI-TWO) DO 165 I=2,M-1 ABTI=AB+TWO*I 165 BJ(I)=I*(ALPHA+I)/(ABTI-ONE)*(BETA+I)/ 1 (ABTI**2)*(AB+I)/(ABTI+ONE) BJ(M)=ZERO 170 DO 180 I=1,M BJ(I)= SQRT(BJ(I)) 180 CONTINUE RETURN END SUBROUTINE WM(W,M,KIND,ALPHA,BETA,IER) c*********************************************************************72 C C ROUTINE TO EVALUATE THE FIRST M MOMENTS OF CLASSICAL WEIGHT C FUNCTIONS C INPUT AND OUTPUT VARIABLES - C O I I I I O C SUBROUTINE WM(W,M,KIND,ALPHA,BETA,IER) C C FUNCTIONS AND SUBROUTINES REFERENCED - DGAMMA SQRT PARCHK DOUBLE PRECISION ALPHA,ALS,BETA,SUM,TMPA,TMPB,TRM,W,DGAMMA DOUBLE PRECISION ZERO,HALF,ONE,TWO DOUBLE PRECISION THREE,FOUR DOUBLE PRECISION PI INTEGER I,IER,JA,JB,K,KIND,M DIMENSION W(M) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) PARAMETER (THREE=3.0D0,FOUR=4.0D0) PARAMETER (PI=3.14159 26535 89793 23846 26433 83279 50 D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) CS PARAMETER (THREE=3.0E0,FOUR=4.0E0) CS PARAMETER (PI=3.14159 26535 89793 23846 26433 83279 50 E0) C CALL PARCHK(KIND,M,ALPHA,BETA,IER) IF(IER.NE.0) RETURN DO 10 K=2,M,2 10 W(K)=ZERO C LEG,CHEB,GEG, JAC,LAG,HERM,EXP,RAT GO TO( 30, 60 , 80, 100,160, 180,20 ,200),KIND 20 ALS=ALPHA GOTO 40 30 ALS=ZERO 40 DO 50 K=1,M,2 50 W(K)=TWO/(K+ALS) RETURN 60 W(1)=PI DO 70 K=3,M,2 70 W(K)=W(K-2)*(K-TWO)/(K-ONE) RETURN 80 W(1)=SQRT(PI)*DGAMMA(ALPHA+ONE)/DGAMMA(ALPHA+THREE/TWO) DO 90 K=3,M,2 90 W(K)=W(K-2)*(K-TWO)/(TWO*ALPHA+K) RETURN 100 ALS=ALPHA+BETA+ONE W(1)=TWO**ALS*DGAMMA(ALPHA+ONE)/DGAMMA(ALS+ONE) 1 *DGAMMA(BETA+ONE) DO 150 K=2,M SUM=ZERO TRM=ONE DO 130 I=0,(K-2)/2 TMPA=TRM DO 110 JA=1,2*I 110 TMPA=TMPA*(ALPHA+JA)/(ALS+JA) DO 120 JB=1,K-2*I-1 120 TMPA=TMPA*(BETA+JB)/(ALS+2*I+JB) TMPA=TMPA/(2*I+ONE)* 1 (2*I*(BETA+ALPHA)+BETA-(K-ONE)*ALPHA)/(BETA+K-2*I-ONE) SUM=SUM+TMPA TRM=TRM*(K-2*I-ONE)/(2*I+ONE)*(K-2*I-TWO)/(2*I+TWO) 130 CONTINUE IF(MOD(K,2).EQ.0) GOTO 150 TMPB=ONE DO 140 I=1,K-1 140 TMPB=TMPB*(ALPHA+I)/(ALS+I) SUM=SUM+TMPB 150 W(K)=SUM*W(1) RETURN 160 W(1)=DGAMMA(ALPHA+ONE) DO 170 K=2,M 170 W(K)=(ALPHA+K-ONE)*W(K-1) RETURN 180 W(1)=DGAMMA((ALPHA+ONE)/TWO) DO 190 K=3,M,2 190 W(K)=W(K-2)*(ALPHA+K-TWO)/TWO RETURN 200 W(1)=DGAMMA(ALPHA+ONE)*DGAMMA(-ALPHA-BETA-ONE)/DGAMMA(-BETA) DO 210 K=2,M 210 W(K)=-W(K-1)*(ALPHA+K-ONE)/(ALPHA+BETA+K) RETURN END SUBROUTINE PARCHK(KIND,M,ALPHA,BETA,IER) c*********************************************************************72 c C ROUTINE TO CHECK RANGES OF PARAMETERS ALPHA, BETA FOR CLASSICAL C WEIGHT FUNCTIONS. M IS THE ORDER OF THE JACOBI MATRIX C REQUIRED AND IS CONSTRAINED BY ALPHA AND BETA FOR THE RATIONAL C WEIGHT FUNCTION (SEE BELOW). M CAN BE REPLACED BY A DUMMY FOR C OTHER WEIGHT FUNCTIONS. C INPUT AND OUTPUT VARIABLES - C I I I I O C SUBROUTINE PARCHK(KIND,M,ALPHA,BETA,IER) C C FUNCTIONS AND SUBROUTINES REFERENCED - WM DOUBLE PRECISION ALPHA,BETA,TMP DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER IER,KIND,M PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C CONSTRAINTS ON ALPHA,BETA:- C 1 ALPHA.GT.-1 C 2 FOR KIND.LT.8 NEED BETA.GT.-1 C 3 FOR KIND.EQ.8 NEED BETA.LT.(ALPHA+BETA+2*M).LT.0 TO C COMPUTE M ELEMENTS OF THE JACOBI MATRIX. C INPUT: C KIND...1-8 FOR CLASSICAL WEIGHT FUNCTION, 0 FOR UNKNOWN) C ALPHA,BETA...AS IN CLASS C M...ORDER OF HIGHEST MOMENT TO BE CALCULATED C OUTPUT: C IER...ERROR INDICATOR - CODED AS FOLLOWS C 1...ALPHA.LE.-1 C 2...BETA.LE.-1 C 3...ALPHA,BETA COMBINATION WRONG FOR RATIONAL WEIGHT C FUNCTION C 4...KIND=0. PARAMETERS CANNOT BE CHECKED AND JACOBI MATRIX C IS NOT OF CLASSICAL TYPE IER=0 IF(KIND.LE.0)IER=4 C CHECK GEGENBAUER,JACOBI,LAGUERRE,HERMITE, EXPONENTIAL IF(KIND.GE.3.AND.(ALPHA.LE.-ONE))IER=1 C C CHECK BETA FOR JACOBI IF(KIND.EQ.4.AND.BETA.LE.-ONE)IER=2 C C CHECK RANGE FOR RATIONAL IF(KIND.LT.8) RETURN TMP=ALPHA+BETA+M+ONE IF(TMP.GE.ZERO.OR.TMP.LE.BETA)IER=3 RETURN END SUBROUTINE CHKQFS(T,WTS,MLT,NT,NWTS,NDX,KEY,W,MOP,MEX, 1 KIND,ALPHA,BETA,LO,E,ER,QM,IER) c*********************************************************************72 C C ROUTINE TO CHECK THE POLYNOMIAL ACCURACY OF A QUADRATURE FORMULA. C IT WILL OPTIONALLY PRINT WEIGHTS, AND RESULTS OF A MOMENTS TEST. C C INPUT AND OUTPUT VARIABLES - C I I I I I I I * I I C SUBROUTINE CHKQFS(T,WTS,MLT,NT,NWTS,NDX,KEY,W,MOP,MEX, C 1 KIND,ALPHA,BETA,LO,E,ER,QM,IER) C I I I I O O O O C C T...ARRAY OF DISTINCT KNOTS C W...MOMENTS ARRAY OF LENGTH MEX C MOP...THE EXPECTED ORDER OF PRECISION OF THE QF C MEX...THE TOTAL NUMBER (.GT.1) OF MOMENTS REQUIRED TO BE TESTED C SET MEX=1 AND LO.LT.0 FOR NO MOMENTS CHECK C KIND...KIND OF CLASSICAL FORMULA. C KIND=0 MEANS UNKNOWN WEIGHT FUNCTION. C THE FIRST MEX MOMENTS MUST BE SET UP C IN ARRAY W BY THE USER FOR THIS CASE. C LO...PRINTING (IF ANY) IS DONE ON UNIT ABS(LO). LO IS CODED C AS FOLLOWS:- C LO.GT.0 MEANS PRINT WEIGHTS AND MOMENT TESTS C LO.EQ.0 MEANS PRINT NOTHING. COMPUTE MOMENT TEST C LO.LT.0 MEANS PRINT WEIGHTS ONLY. DON'T COMPUTE MOMENT TESTS C E,ER...ABSOLUTE AND RELATIVE ERRORS OF THE QF APPLIED C TO (X-DEL)**N C QM...VALUES OF THE QF APPLIED TO (X-DEL)**N C IER...ERROR INDICATOR. ANY ERROR RETURN COMES FROM WM. C C FUNCTIONS AND SUBROUTINES REFERENCED - WM DOUBLE PRECISION ALPHA,BETA,E,EK,EMN,EMX,EREST,ERN,ERX,ER,PX DOUBLE PRECISION TMP,TMPX,PREC,QM,T,W,WTS DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,J,JL,K,KEY,KIND,KINDP,KJL,L,LO,LU,M,MEX,MLT INTEGER MOP,MX,NDX,NT,NWTS DIMENSION T(NT),WTS(NWTS),MLT(NT),W(MEX),NDX(NT) DIMENSION E(MEX),ER(MEX),QM(MEX) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) COMMON /CTRLR/ PREC(10) CHARACTER*72 TXT1(10) DATA TXT1/ 1' INTERPOLATORY QUADRATURE FORMULA', 2' TYPE INTERVAL WEIGHT FUNCTION NAME', 3' 1 (-1,1) ONE LEGENDRE', 4' 2 (-1,1) (1-X**2)**(-HALF) CHEBYSHEV', 5' 3 (-1,1) (1-X**2)**ALPHA GEGENBAUER', 6' 4 (-1,1) (1-X)**ALPHA*(1+X)**BETA JACOBI', 7' 5 (0,INF) X**ALPHA*EXP(-X) GEN LAGUERRE', 8' 6 (-INF,INF) ABS(X)**ALFA*EXP(-X**2) GEN HERMITE', 9' 7 (-1,1) ABS(X)**ALFA EXPONENTIAL', 1' 8 (0,INF) X**ALFA*(1+X)**BETA RATIONAL'/ LU=ABS(LO) C KIND MAY BE SET TO -1 TO ALLOW PRINTING OF MOMENTS ONLY C THIS FEATURE IS ONLY USED INTERNALLY (BY CHKQF) KINDP=MAX(0,KIND) IF(LO.EQ.0.OR.KIND.EQ.-1) GOTO 40 C IF(KINDP.EQ.0) GOTO 10 WRITE (LU,150)TXT1(1),TXT1(2),TXT1(KINDP+2) IF(KINDP.GE.3) WRITE (LU,160) ALPHA IF(KINDP.EQ.4.OR.KINDP.EQ.8) WRITE (LU,170) BETA 10 IF(KIND.NE.-1) WRITE(LU,180)PREC(1) C WRITE(LU,210) DO30 I=1,NT K=ABS(NDX(I)) IF(K.EQ.0) GOTO 30 WRITE(LU,190)I,T(I),MLT(I),WTS(K) DO 20 J=K+1,K+MLT(I)-1 WRITE(LU,200)WTS(J) 20 CONTINUE 30 CONTINUE 40 IF(LO.LT.0) RETURN IER=0 IF(KINDP.NE.0)CALL WM(W,MEX,KINDP,ALPHA,BETA,IER) IF(IER.NE.0) RETURN 50 DO 60 K=1,MEX QM(K)=ZERO 60 CONTINUE EREST=ZERO DO 90 K=1,NT TMP=ONE L=ABS(NDX(K)) IF(L.EQ.0) GOTO 90 EREST=EREST+ABS(WTS(L)) DO 80 J=1,MEX QM(J)=QM(J)+TMP*WTS(L) TMPX=TMP PX=ONE DO 70 JL=2,MIN(MLT(K),MEX-J+1) KJL=J+JL-1 TMPX=TMPX*(KJL-1) QM(KJL)=QM(KJL)+TMPX*WTS(L+JL-1)/PX IF(KEY.GT.0) GOTO 70 PX=PX*JL 70 CONTINUE TMP=TMP*T(K) 80 CONTINUE 90 CONTINUE DO 100 K=1,MEX E(K)=W(K)-QM(K) ER(K)=E(K)/(ABS(W(K))+ONE) 100 CONTINUE C FOR SOME STRANGE WEIGHT FUNCTIONS W(1) MAY VANISH EREST=EREST/(ABS(W(1))+ONE) C EXIT IF USER DOES NOT WANT PRINTED OUTPUT IF(LO.EQ.0) RETURN EMX=ABS(E(1)) EMN=EMX ERX=ABS(ER(1)) ERN=ERX M=MOP+1 MX=MIN(MOP,MEX) DO 110 K=2,MX EMX=MAX(ABS(E(K)),EMX) EMN=MIN(ABS(E(K)),EMN) ERX=MAX(ABS(ER(K)),ERX) ERN=MIN(ABS(ER(K)),ERN) 110 CONTINUE WRITE(LU,220) MOP,EMN,ERN,EMX,ERX,EREST IF(MEX.LT.M) GOTO 140 120 EK=E(M) DO 130 J=1,MOP EK=EK/J 130 CONTINUE WRITE(LU,230) MOP,E(M),EK 140 WRITE (LU,240) WRITE (LU,250) (J,W(J),QM(J),E(J),ER(J),J=1,MX) WRITE (LU,250) IF(MEX.GE.M)WRITE (LU,250) (J,W(J),QM(J),E(J),ER(J),J=M,MEX) 150 FORMAT(1H1,(8X,A72/)) 160 FORMAT(/24H PARAMETER(S) ALPHA ,F12.5) 170 FORMAT( 24H BETA ,F12.5) 180 FORMAT(/23H MACHINE PRECISION ,D13.1) 190 FORMAT(2(I4,D26.17)) 200 FORMAT(37X,D26.17) 210 FORMAT(/11X,5HKNOTS,15X,10HMULT ,10X,7HWEIGHTS/) 220 FORMAT(//' COMPARISON OF MOMENTS',// 1,' ORDER OF PRECISION',I4,// 2,' ERRORS : ABSOLUTE RELATIVE'/ 3,'---------+-------------------------'/ 4,' MINIMUM :',2D12.3/ 5,' MAXIMUM :',2D12.3// 6,' WEIGHTS RATIO ',D13.3) 230 FORMAT(' ERROR FOR ',I3,'-TH POWER ',D13.3/ 1 ,' ERROR CONSTANT ',D13.3,/) 240 FORMAT (/12H MOMENTS: /12X,4HTRUE ,12X, 1 9HFROM Q.F. ,9X,5HERROR ,5X,8HRELATIVE /) 250 FORMAT (I4,2D19.10,2D12.3) RETURN END SUBROUTINE CHKQF(T,WTS,MLT,NT,NWTS,NDX,KEY,WF,MOP,MEX,KIND, 1 ALPHA,BETA,LO,E,ER,QM,NWF,A,B,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE AND PRINT THE MOMENTS OF A QF FOR C A CLASICAL WEIGHT FUNCTION WITH ANY VALID A,B C NO CHECK CAN BE MADE FOR NON-CLASSICAL WEIGHT FUNCTIONS C C INPUT AND OUTPUT VARIABLES - C I I I I I I I O I I I C SUBROUTINE CHKQF(T,WTS,MLT,NT,NWTS,NDX,KEY,WF,MOP,MEX,KIND, C 1 ALPHA,BETA,LO,E,ER,QM,NWF,A,B,IER) C I I I O O O I I I O C C NWF...SIZE OF WORKFIELD ARRAY. MUST BE .GE. MEX+NT C MOP...THE EXPECTED ORDER OF PRECISION OF THE QF. C MEX...THE TOTAL NUMBER (.GT.1) OF MOMENTS REQUIRED TO BE TESTED C SET MEX=1 AND LO.LT.0 FOR NO MOMENTS CHECK C LO...PRINTING (IF ANY) IS DONE ON UNIT ABS(LO). LO IS CODED C AS FOLLOWS:- C LO.GT.0 MEANS PRINT WEIGHTS AND MOMENT TESTS C LO.EQ.0 MEANS PRINT NOTHING. COMPUTE MOMENT TEST C LO.LT.0 MEANS PRINT WEIGHTS ONLY. DON'T COMPUTE MOMENT TESTS C E,ER...ABSOLUTE AND RELATIVE ERRORS OF THE QF APPLIED C TO (X-DEL)**N C QM...VALUES OF THE QF APPLIED TO (X-DEL)**N C IER...ERROR INDICATOR. ANY ERROR RETURN COMES FROM WM. C C IER CODES - C 1-4 ERROR RETURN FROM PARCHK: ALPHA, BETA WRONG C SEE ROUTINE PARCHECK C 6 ZERO LENGTH INTERVAL (KIND=1,2,3,4,7) C 7 B.LE.0 FOR (KIND=5,6) C 8 A+B.LE.0 (KIND=8) C C FUNCTIONS AND SUBROUTINES REFERENCED - CHKQFS PARCHK SCMM DOUBLE PRECISION A,ALPHA,B,BETA,E,ER,PREC,QM,T,TMP,WF,WTS DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,IZERO,KEY,KIND,LEX,LO,LU,MEX,MOP,NA,NEG,NT,NWF INTEGER NWTS,MLT,NDX DIMENSION T(NT),MLT(NT),WTS(NWTS),NDX(NT),WF(NWF) DIMENSION E(MEX),ER(MEX),QM(MEX) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) COMMON /CTRLR/ PREC(10) CHARACTER*72 TXT2(10) DATA TXT2/ 1' INTERPOLATORY QUADRATURE FORMULA', 2' TYPE INTERVAL WEIGHT FUNCTION NAME', 3' 1 (A,B) ONE LEGENDRE', 4' 2 (A,B) ((B-X)*(X-A))**(-HALF) CHEBYSHEV', 5' 3 (A,B) ((B-X)*(X-A))**ALPHA GEGENBAUER', 6' 4 (A,B) (B-X)**ALPHA*(X-A)**BETA JACOBI', 7' 5 (A,INF) (X-A)**ALPHA*EXP(-B*(X-A)) GEN LAGUERRE', 8' 6 (-INF,INF) ABS(X-A)**ALFA*EXP(-B*(X-A)**2) GEN HERMITE', 9' 7 (A,B) ABS(X-(A+B)/TWO)**ALFA EXPONENTIAL', 1' 8 (A,INF) (X-A)**ALFA*(B+X)**BETA RATIONAL'/ C CALL PARCHK(KIND,MEX,ALPHA,BETA,IER) IF(IER.NE.0) RETURN IF(LO.EQ.0) GOTO 10 C PRINT WEIGHTS IZERO=0 LU=ABS(LO) WRITE (LU,50)TXT2(1),TXT2(2),TXT2(KIND+2) WRITE (LU,80) A WRITE (LU,90) B IF (KIND.GE.3) WRITE (LU,60) ALPHA IF (KIND.EQ.4.OR.KIND.EQ.8) WRITE (LU,70) BETA CALL CHKQFS(T,WTS,MLT,NT,NWTS,NDX,KEY,WF,MOP,MEX,IZERO, 1 ALPHA,BETA,-LU,E,ER,QM,IER) IF (IER.NE.0.OR.LO.LT.0) RETURN 10 LEX=MEX+NT IF(NWF.GE.LEX) GOTO 20 IER=-LEX RETURN 20 CONTINUE CALL SCMM(WF,MEX,KIND,ALPHA,BETA,A,B,IER) IF (IER.NE.0) RETURN 30 NA=MEX+1 TMP=(B+A)/TWO IF(KIND.EQ.5.OR.KIND.EQ.6.OR.KIND.EQ.8)TMP=A DO 40 I=1,NT WF(NA+I-1)=T(I)-TMP 40 CONTINUE NEG=-1 C CHECK MOMENTS CALL CHKQFS(WF(NA),WTS,MLT,NT,NWTS,NDX,KEY,WF,MOP,MEX,NEG, 1 ALPHA,BETA,LO,E,ER,QM,IER) RETURN 50 FORMAT(1H1,(8X,A72/)) 60 FORMAT ( ' ALPHA ',F12.5) 70 FORMAT ( ' BETA ',F12.5) 80 FORMAT ( ' PARAMETERS A ',F12.5) 90 FORMAT ( ' B ',F12.5) 100 FORMAT (/23H MACHINE PRECISION ,D13.1) END SUBROUTINE SCT(NT,T,ST,KIND,A,B,IER) c*********************************************************************72 C C ROUTINE TO SCALE DISTINCT KNOTS FOR ANY VALID A,B TO THOSE FOR C THE DEFAULT VALUES OF A,B. ARRAYS T AND ST MAY COINCIDE. C ALL KNOTS IN THE T ARRAY ARE SCALED AND ARE OUTPUT IN ST. C C INPUT AND OUTPUT VARIABLES - C I I O I I I O C SUBROUTINE SCT(NT,T,ST,KIND,A,B,IER) C C FUNCTIONS AND SUBROUTINES REFERENCED - MACHEP DOUBLE PRECISION A,B,BMA,SHFT,SLP,TMP,ST,T DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,KIND,NT DIMENSION T(NT),ST(NT) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C IER=0 IF(KIND.LE.0.OR.KIND.GT.8) THEN IER=4 RETURN ENDIF C LEG,CHEB,GEG,JAC,LAG,HERM,EXP,RAT 10 GO TO( 20, 20, 20, 20, 30, 40, 20, 50),KIND 20 CONTINUE CALL MACHEP(TMP) BMA=B-A IF(BMA.GT.TMP) THEN SLP=TWO/BMA SHFT=-(A+B)/BMA ELSE IER=6 RETURN ENDIF GOTO 60 30 CONTINUE IF(B.LT.ZERO) THEN IER=7 RETURN ENDIF SLP=B SHFT=-A*B GOTO 60 40 CONTINUE IF(B.LT.ZERO) THEN IER=7 RETURN ENDIF SLP=SQRT(B) SHFT=-A*SLP GOTO 60 50 CONTINUE SLP=ONE/(A+B) IF(SLP.LE.ZERO) THEN IER=8 RETURN ENDIF SHFT=-A*SLP 60 CONTINUE DO 70 I=1,NT ST(I)=SHFT+SLP*T(I) 70 CONTINUE RETURN END SUBROUTINE SCQF(NT,T,MLT,WTS,NWTS,NDX,SWTS,ST, 1 KIND,ALPHA,BETA,A,B,IER) c*********************************************************************72 c C ROUTINE TO SCALE WEIGHTS AND KNOTS FOR CLASSICAL WEIGHT FUNCTION C WITH DEFAULT VALUES FOR A AND B TO THOSE FOR ANY VALID A,B C C INPUT AND OUTPUT VARIABLES - C I I I I I I O O C SUBROUTINE SCQF(NT,T,MLT,WTS,NWTS,NDX,SWTS,ST, C 1 KIND,ALPHA,BETA,A,B,IER) C I I I I I O C C THE ARRAYS WTS AND SWTS MAY COINCIDE C THE ARRAYS T AND ST MAY COINCIDE C IER CODES - C 1-4 ERROR RETURN FROM PARCHK: ALPHA, BETA WRONG C SEE ROUTINE PARCHECK C 6 ZERO LENGTH INTERVAL (KIND=1,2,3,4,7) C 7 B.LE.0 FOR (KIND=5,6) C 8 A+B.LE.0 (KIND=8) C C C FUNCTIONS AND SUBROUTINES REFERENCED - MACHEP SQRT PARCHK DOUBLE PRECISION A,AL,ALPHA,B,BE,BETA,P,SHFT,SLP,TEMP DOUBLE PRECISION TMP,ST,SWTS,T,WTS DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,K,KIND,L,NT,NWTS,MLT,NDX DIMENSION T(NT),MLT(NT),WTS(NWTS),NDX(NT),SWTS(NWTS),ST(NT) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C CALL MACHEP(TEMP) CALL PARCHK(KIND,1,ALPHA,BETA,IER) IF(IER.NE.0) RETURN C C LEG,CHEB,GEG,JAC,LAG,HERM,EXP,RAT 10 GO TO( 20, 30, 40, 50, 90, 110, 60,130),KIND 20 AL=ZERO BE=ZERO GOTO 70 30 AL=-HALF BE=-HALF GOTO 70 40 AL=ALPHA BE=ALPHA GOTO 70 50 AL=ALPHA BE=BETA GOTO 70 60 AL=ALPHA BE=ZERO 70 IF((B-A).GT.TEMP) GOTO 80 IER=6 RETURN 80 SHFT=(A + B)/TWO SLP=(B - A)/TWO GOTO 150 90 IF(B.GT.ZERO) GOTO 100 IER=7 RETURN 100 SHFT=A SLP=ONE/B AL=ALPHA BE=ZERO GOTO 150 110 IF(B.GT.ZERO) GOTO 120 IER=7 RETURN 120 SHFT=A SLP=ONE/SQRT(B) AL=ALPHA BE=ZERO GOTO 150 130 IF(A+B.GT.ZERO) GOTO 140 IER=8 RETURN 140 SHFT=A SLP=A + B AL=ALPHA BE=BETA 150 CONTINUE P=SLP**(AL+BE+ONE) DO 170 K=1,NT ST(K)=SHFT + SLP*T(K) L=ABS(NDX(K)) IF(L.EQ.0) GOTO 170 TMP=P DO 160 I=L,L+MLT(K)-1 SWTS(I)=WTS(I)*TMP 160 TMP=TMP*SLP 170 CONTINUE RETURN END SUBROUTINE SCMM(W,M,KIND,ALPHA,BETA,A,B,IER) c*********************************************************************72 c C ROUTINE TO COMPUTE MOMENTS OF CLASSICAL WEIGHT FUNCTION WITH C DEFAULT VALUES FOR A,B AND SCALE THEM TO THOSE FOR ANY VALID A,B C INPUT AND OUTPUT VARIABLES - C O I I I I I I O C SUBROUTINE SCMM(W,M,KIND,ALPHA,BETA,A,B,IER) C C MOMENTS ARE OUTPUT TO W C FUNCTIONS AND SUBROUTINES REFERENCED - MACHEP SQRT WM DOUBLE PRECISION A,AL,ALPHA,B,BE,BETA,P,Q,TEMP,TMP,W DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,KIND,M DIMENSION W(M) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C CALL MACHEP(TEMP) C C LEG,CHEB,GEG,JAC,LAG,HERM,EXP, RAT 10 GO TO( 20, 30, 40, 50, 90, 110, 60, 130),KIND 20 AL=ZERO BE=ZERO GOTO 70 30 AL=-HALF BE=-HALF GOTO 70 40 AL=ALPHA BE=ALPHA GOTO 70 50 AL=ALPHA BE=BETA GOTO 70 60 AL=ALPHA BE=ZERO 70 IF((B-A).GT.TEMP) GOTO 80 IER=6 RETURN 80 Q=(B-A)/TWO P=Q**(AL+BE+ONE) GOTO 150 90 IF(B.GT.ZERO) GOTO 100 IER=7 RETURN 100 Q=ONE/B P=Q**(ALPHA+ONE) GOTO 150 110 IF(B.GT.ZERO) GOTO 120 IER=7 RETURN 120 Q=ONE/SQRT(B) P=Q**(ALPHA+ONE) GOTO 150 130 IF(A+B.GT.ZERO) GOTO 140 IER=8 RETURN 140 CONTINUE Q=A+B P=Q**(ALPHA+BETA+ONE) 150 CALL WM(W,M,KIND,ALPHA,BETA,IER) IF(IER.EQ.0) GOTO 160 RETURN 160 TMP=P DO 170 I=1,M W(I)=W(I)*TMP TMP=TMP*Q 170 CONTINUE RETURN END SUBROUTINE WTFN(T,W,NT,KIND,ALPHA,BETA,IER) c*********************************************************************72 C C ROUTINE TO EVALUATE THE CLASSICAL WEIGHT FUNCTIONS AT THE POINTS C GIVEN IN ARRAY T. THE INPUT, T, AND OUTPUT, W, ARRAYS MAY BE THE C SAME. C INPUT AND OUTPUT VARIABLES - C I O I I I I O C SUBROUTINE WTFN(T,W,NT,KIND,ALPHA,BETA,IER) C C *******WARNING******* C NO CHECK IS MADE C (1) THAT THE WEIGHT FUNCTION IS DEFINED FOR THE POINTS IN T. C (2) THAT THE POINTS ARE IN THE APPROPRIATE INTERVAL. C HOWEVER THE PARAMETERS ALPHA AND BETA ARE CHECKED FOR THE C CLASSICAL WEIGHT FUNCTIONS. C C FUNCTIONS AND SUBROUTINES REFERENCED - EXP SQRT PARCHK DOUBLE PRECISION ALPHA,BETA,T,W DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER IER,K,KIND,NT DIMENSION W(NT),T(NT) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C CALL PARCHK(KIND,1,ALPHA,BETA,IER) IF(IER.NE.0) RETURN C LEG,CHEB,GEG,JAC,LAG,HERM,EXP,RAT GO TO( 30, 50 , 70, 90,130, 160,10 ,190),KIND 10 IF (ALPHA.EQ.ZERO) GOTO 30 DO 20 K=1,NT 20 W(K)=ABS(T(K))**ALPHA RETURN 30 DO 40 K=1,NT 40 W(K)=ONE RETURN 50 DO 60 K=1,NT 60 W(K)=ONE/SQRT((ONE-T(K))*(ONE+T(K))) RETURN 70 IF (ALPHA.EQ.ZERO) GOTO 30 DO 80 K=1,NT 80 W(K)=((ONE-T(K))*(ONE+T(K)))**ALPHA RETURN 90 DO 100 K=1,NT 100 W(K)=ONE IF (ALPHA.NE.ZERO) THEN DO 110 K=1,NT 110 W(K)=W(K)*(ONE-T(K))**ALPHA END IF IF (BETA.NE.ZERO) THEN DO 120 K=1,NT 120 W(K)=W(K)*(ONE+T(K))**BETA END IF RETURN 130 DO 140 K=1,NT 140 W(K)=EXP(-T(K)) IF (ALPHA.NE.ZERO) THEN DO 150 K=1,NT 150 W(K)=W(K)*T(K)**ALPHA END IF RETURN 160 DO 170 K=1,NT 170 W(K)=EXP(-T(K)**2) IF (ALPHA.NE.ZERO) THEN DO 180 K=1,NT 180 W(K)=W(K)*ABS(T(K))**ALPHA END IF RETURN 190 DO 200 K=1,NT 200 W(K)=ONE IF (ALPHA.NE.ZERO) THEN DO 210 K=1,NT 210 W(K)=W(K)*T(K)**ALPHA END IF IF (BETA.NE.ZERO) THEN DO 220 K=1,NT 220 W(K)=W(K)*(ONE+T(K))**BETA END IF RETURN END SUBROUTINE IMTQLX(N,D,E,Z,IER) c*********************************************************************72 c C THIS ROUTINE IS A SLIGHTLY MODIFIED VERSION OF THE EISPACK C ROUTINE TO PERFORM THE IMPLICIT QL ALGORITHM ON A SYMMETRIC C TRIDIAGONAL MATRIX. THE AUTHORS THANK THE AUTHORS OF EISPACK C FOR PERMISSION TO USE THIS ROUTINE. FOR DETAILS SEE C MARTIN AND WILKINSON: THE IMPLICIT QL ALGORITHM, NUMER MATH C 12, 277-383 (1968). IT HAS BEEN MODIFIED TO PRODUCE THE C T C PRODUCT Q Z, WHERE Z IS AN INPUT VECTOR AND Q IS THE C ORTHOGONAL MATRIX DIAGONALIZING THE INPUT MATRIX. THE CHANGES C CONSIST (ESSENTIALLY) OF APPLYING THE ORTHOGONAL TRANSFORMATIONS C DIRECTLY TO Z AS THEY ARE GENERATED. SEE REFERENCES TO Z NEAR C STATEMENT 60. C C FUNCTIONS AND SUBROUTINES REFERENCED - SIGN SQRT DOUBLE PRECISION B,C,F,G,P,R,S,D,E,PREC,Z DOUBLE PRECISION ZERO,HALF,ONE,TWO INTEGER I,IER,II,J,K,L,M,MML,N INTEGER ITN DIMENSION D(N),E(N),Z(N) PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) PARAMETER (ITN=30) COMMON/CTRLR/PREC(10) IER=0 IF(N.EQ.1) GOTO 110 E(N)=ZERO DO 70 L=1,N J=0 10 DO 20 M=L,N IF(M.EQ.N) GOTO 30 IF( ABS(E(M)).LE.PREC(1)*( ABS(D(M))+ ABS(D(M+1)))) GOTO 30 20 CONTINUE 30 P=D(L) IF(M.EQ.L) GOTO 70 IF(J.EQ.ITN) GOTO 100 J=J+1 G=(D(L+1)-P)/(TWO*E(L)) R= SQRT(G*G+ONE) G=D(M)-P+E(L)/(G+ SIGN(R,G)) S=ONE C=ONE P=ZERO MML=M-L DO 60 II=1,MML I=M-II F=S*E(I) B=C*E(I) IF( ABS(F).LT. ABS(G)) GOTO 40 C=G/F R= SQRT(C*C+ONE) E(I+1)=F*R S=ONE/R C=C*S GOTO 50 40 S=F/G R= SQRT(S*S+ONE) E(I+1)=G*R C=ONE/R S=S*C 50 G=D(I+1)-P R=(D(I)-G)*S+TWO*C*B P=S*R D(I+1)=G+P G=C*R-B F=Z(I+1) Z(I+1)=S*Z(I)+C*F 60 Z(I)=C*Z(I)-S*F D(L)=D(L)-P E(L)=G E(M)=ZERO GOTO 10 70 CONTINUE DO 90 II=2,N I=II-1 K=I P=D(I) DO 80 J=II,N IF(D(J).GE.P) GOTO 80 K=J P=D(J) 80 CONTINUE IF(K.EQ.I) GOTO 90 D(K)=D(I) D(I)=P P=Z(I) Z(I)=Z(K) Z(K)=P 90 CONTINUE GOTO 110 100 IER=L write ( *, * ) ' ' write ( *, * ) 'IMTQLX - Fatal error!' write ( *, * ) ' Iteration limit exceeded.' 110 RETURN END SUBROUTINE MACHEP (X) c*********************************************************************72 c C ROUTINE TO COMPUTE MACHINE EPSILON, C DEFINED AS THE SMALLEST FLOATING POINT NUMBER SUCH THAT C FLOAT(1.0+EPS) > 1.0. OUTPUT GOES TO X DOUBLE PRECISION X,Y DOUBLE PRECISION ZERO,HALF,ONE,TWO PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) CS PARAMETER (ZERO=0.0E0,HALF=0.5E0,ONE=1.0E0,TWO=2.0E0) C Y=ONE 10 IF (Y+ONE.NE.ONE) THEN X=Y Y=Y/TWO GOTO 10 END IF RETURN END DOUBLE PRECISION FUNCTION DGAMMA(X) DGA 10 CS REAL FUNCTION GAMMA(X) DGA 20 c*********************************************************************72 C DGA 40 C THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. DGA 50 C COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN W. J. CODY, DGA 60 C 'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL FUNCTIONS', DGA 70 C LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS DUNDEE, DGA 80 C 1975, G. A. WATSON (ED.), SPRINGER VERLAG, BERLIN, 1976. THE DGA 90 C PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA DGA 100 C FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS DGA 110 C FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. DGA 120 C THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM HART, ET. AL., DGA 130 C COMPUTER APPROXIMATIONS, WILEY AND SONS, NEW YORK, 1968. LOWER DGA 140 C ORDER APPROXIMATIONS CAN BE SUBSTITUTED FOR THESE ON MACHINES DGA 150 C WITH LESS PRECISE ARITHMETIC. DGA 160 C DGA 170 C DGA 180 C******************************************************************* DGA 190 C******************************************************************* DGA 200 C DGA 210 C EXPLANATION OF MACHINE-DEPENDENT CONSTANTS DGA 220 C DGA 230 C EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT DGA 240 C 1.0 + EPS .GT. 1.0 DGA 250 C XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE DGA 260 C IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION DGA 270 C GAMMA(XBIG) = XINF. DGA 280 C XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT DGA 290 C 1/XMININ IS MACHINE REPRESENTABLE. DGA 300 C XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER. DGA 310 C DGA 320 C APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: DGA 330 C DGA 340 C IBM/195 CDC/7600 UNIVAC/1108 VAX 11/780 (UNIX) DGA 350 C (D.P.) (S.P.,RNDG) (D.P.) (S.P.) (D.P.) DGA 360 C DGA 370 C EPS 2.221D-16 3.553E-15 1.735D-18 5.961E-08 1.388D-16 DGA 380 C XBIG 57.574 177.802 171.489 34.844 34.844 DGA 390 C XMININ 1.382D-76 3.132E-294 1.113D-308 5.883E-39 5.883D-39 DGA 400 C XINF 7.23D+75 1.26E+322 8.98D+307 1.70E+38 1.70D+38 DGA 410 C DGA 420 C******************************************************************* DGA 430 C******************************************************************* DGA 440 C DGA 450 C DGA 460 C ERROR RETURNS DGA 470 C DGA 480 C THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR DGA 490 C WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED DGA 500 C TO BE FREE OF UNDERFLOW AND OVERFLOW. DGA 510 C DGA 520 C DGA 530 C DGA 540 C OTHER SUBPROGRAMS REQUIRED (SINGLE PRECISION VERSION) DGA 550 C DGA 560 C ALOG,EXP,FLOAT,IFIX,SIN DGA 570 C DGA 580 C OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) DGA 590 C DGA 600 C DBLE,DEXP,DLOG,DSIN,FLOAT,IFIX,SNGL DGA 610 C DGA 620 C DGA 630 C DGA 640 C AUTHOR: W. J. CODY DGA 650 C APPLIED MATHEMATICS DIVISION DGA 660 C ARGONNE NATIONAL LABORATORY DGA 670 C ARGONNE, IL 60439 DGA 680 C DGA 690 C LATEST MODIFICATION: MAY 18, 1982 DGA 700 C DGA 710 C---------------------------------------------------------------------- DGA 720 CS REAL C,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI, DGA 730 CS 1 SUM,TWELVE,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DGA 740 DOUBLE PRECISION C, EPS, FACT, HALF, ONE, P, PI, Q, RES, SQRTPI, DGA 750 * SUM, TWELVE, X, XBIG, XDEN, XINF, XMININ, XNUM, Y, Y1, YSQ, Z, DGA 760 * ZERO DGA 770 INTEGER I, J, N DGA 780 LOGICAL PARITY DGA 790 DIMENSION C(7), P(8), Q(8) DGA 800 C---------------------------------------------------------------------- DGA 810 C MATHEMATICAL CONSTANTS DGA 820 C---------------------------------------------------------------------- DGA 830 CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/ DGA 840 CS DATA SQRTPI/0.9189385332046727417803297E0/ DGA 850 CS DATA PI/3.1415926535897932384626434E0/ DGA 860 DATA ONE, HALF, TWELVE, ZERO /1.0D0,0.5D0,12.0D0,0.0D0/ DGA 870 DATA SQRTPI /0.9189385332046727417803297D0/ DGA 880 DATA PI /3.1415926535897932384626434D0/ DGA 890 C---------------------------------------------------------------------- DGA 900 C MACHINE DEPENDENT PARAMETERS (VAX 11/780 DP) DGA 910 C---------------------------------------------------------------------- DGA 920 CS DATA XBIG,XMININ,EPS/34.844E0,5.883E-39,5.961E-08/ DGA 930 CS DATA XINF/1.7014E38/ DGA 940 DATA XBIG, XMININ, EPS /34.844D0,5.883D-39,1.388D-17/ DGA 950 DATA XINF /1.7014D38/ DGA 960 C---------------------------------------------------------------------- DGA 970 C NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX DGA 980 C APPROXIMATION OVER (1,2). DGA 990 C---------------------------------------------------------------------- DGA 1000 CS DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1,DGA 1010 CS 1 -3.79804256470945635097577E+2,6.29331155312818442661052E+2,DGA 1020 CS 2 8.66966202790413211295064E+2,-3.14512729688483675254357E+4,DGA 1030 CS 3 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/DGA 1040 CS DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2,DGA 1050 CS 1 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3,DGA 1060 CS 2 2.25381184209801510330112E+4,4.75584627752788110767815E+3,DGA 1070 CS 3 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/DGA 1080 DATA P /-1.71618513886549492533811D+0,2.47656508055759199108314D+1DGA 1090 * ,-3.79804256470945635097577D+2,6.29331155312818442661052D+2, DGA 1100 * 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, DGA 1110 * -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DGA 1120 DATA Q /-3.08402300119738975254353D+1,3.15350626979604161529144D+2DGA 1130 * ,-1.01515636749021914166146D+3,-3.10777167157231109440444D+3, DGA 1140 * 2.25381184209801510330112D+4,4.75584627752788110767815D+3, DGA 1150 * -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ DGA 1160 C---------------------------------------------------------------------- DGA 1170 C COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). DGA 1180 C---------------------------------------------------------------------- DGA 1190 CS DATA C/-1.910444077728E-03,8.4171387781295E-04, DGA 1200 CS 1 -5.952379913043012E-04,7.93650793500350248E-04, DGA 1210 CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, DGA 1220 CS 3 5.7083835261E-03/ DGA 1230 DATA C /-1.910444077728D-03,8.4171387781295D-04, DGA 1240 * -5.952379913043012D-04,7.93650793500350248D-04, DGA 1250 * -2.777777777777681622553D-03,8.333333333333333331554247D-02, DGA 1260 * 5.7083835261D-03/ DGA 1270 C---------------------------------------------------------------------- DGA 1280 PARITY = .FALSE. DGA 1290 FACT = ONE DGA 1300 N = 0 DGA 1310 Y = X DGA 1320 IF (Y.GT.ZERO) GO TO 10 DGA 1330 C---------------------------------------------------------------------- DGA 1340 C ARGUMENT IS NEGATIVE DGA 1350 C---------------------------------------------------------------------- DGA 1360 Y = -X DGA 1370 CS J = IFIX(Y) DGA 1380 J = IFIX(SNGL(Y)) DGA 1390 CS RES = Y - FLOAT(J) DGA 1400 RES = Y - DBLE(FLOAT(J)) DGA 1410 IF (RES.EQ.ZERO) GO TO 100 DGA 1420 IF (J.NE.(J/2)*2) PARITY = .TRUE. DGA 1430 CS FACT = -PI / SIN(PI*RES) DGA 1440 FACT = -PI/DSIN(PI*RES) DGA 1450 Y = Y + ONE DGA 1460 C---------------------------------------------------------------------- DGA 1470 C ARGUMENT IS POSITIVE DGA 1480 C---------------------------------------------------------------------- DGA 1490 10 IF (Y.LT.EPS) GO TO 90 DGA 1500 IF (Y.GE.TWELVE) GO TO 70 DGA 1510 Y1 = Y DGA 1520 IF (Y.GE.ONE) GO TO 20 DGA 1530 C---------------------------------------------------------------------- DGA 1540 C 0.0 .LT. ARGUMENT .LT. 1.0 DGA 1550 C---------------------------------------------------------------------- DGA 1560 Z = Y DGA 1570 Y = Y + ONE DGA 1580 GO TO 30 DGA 1590 C---------------------------------------------------------------------- DGA 1600 C 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY DGA 1610 C---------------------------------------------------------------------- DGA 1620 CS 20 N = IFIX(Y) - 1 DGA 1630 20 N = IFIX(SNGL(Y)) - 1 DGA 1640 CS Y = Y - FLOAT(N) DGA 1650 Y = Y - DBLE(FLOAT(N)) DGA 1660 Z = Y - ONE DGA 1670 C---------------------------------------------------------------------- DGA 1680 C EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 DGA 1690 C---------------------------------------------------------------------- DGA 1700 30 XNUM = ZERO DGA 1710 XDEN = ONE DGA 1720 DO 40 I=1,8 DGA 1730 XNUM = (XNUM+P(I))*Z DGA 1740 XDEN = XDEN*Z + Q(I) DGA 1750 40 CONTINUE DGA 1760 RES = XNUM/XDEN + ONE DGA 1770 IF (Y.EQ.Y1) GO TO 110 DGA 1780 IF (Y1.GT.Y) GO TO 50 DGA 1790 C---------------------------------------------------------------------- DGA 1800 C ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 DGA 1810 C---------------------------------------------------------------------- DGA 1820 RES = RES/Y1 DGA 1830 GO TO 110 DGA 1840 C---------------------------------------------------------------------- DGA 1850 C ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 DGA 1860 C---------------------------------------------------------------------- DGA 1870 50 DO 60 I=1,N DGA 1880 RES = RES*Y DGA 1890 Y = Y + ONE DGA 1900 60 CONTINUE DGA 1910 GO TO 110 DGA 1920 C---------------------------------------------------------------------- DGA 1930 C EVALUATE FOR ARGUMENT .GE. 12.0, DGA 1940 C---------------------------------------------------------------------- DGA 1950 70 IF (Y.GT.XBIG) GO TO 100 DGA 1960 YSQ = Y*Y DGA 1970 SUM = C(7) DGA 1980 DO 80 I=1,6 DGA 1990 SUM = SUM/YSQ + C(I) DGA 2000 80 CONTINUE DGA 2010 SUM = SUM/Y - Y + SQRTPI DGA 2020 CS SUM = SUM + (Y-HALF)*ALOG(Y) DGA 2030 SUM = SUM + (Y-HALF)*DLOG(Y) DGA 2040 CS RES = EXP(SUM) DGA 2050 RES = DEXP(SUM) DGA 2060 GO TO 110 DGA 2070 C---------------------------------------------------------------------- DGA 2080 C ARGUMENT .LT. EPS DGA 2090 C---------------------------------------------------------------------- DGA 2100 90 IF (Y.LT.XMININ) GO TO 100 DGA 2110 RES = ONE/Y DGA 2120 GO TO 110 DGA 2130 C---------------------------------------------------------------------- DGA 2140 C RETURN FOR SINGULARITIES, EXTREME ARGUMENTS, ETC. DGA 2150 C---------------------------------------------------------------------- DGA 2160 CS100 GAMMA = XINF DGA 2170 100 DGAMMA = XINF DGA 2180 GO TO 120 DGA 2190 C---------------------------------------------------------------------- DGA 2200 C FINAL ADJUSTMENTS AND RETURN DGA 2210 C---------------------------------------------------------------------- DGA 2220 110 IF (PARITY) RES = -RES DGA 2230 IF (FACT.NE.ONE) RES = FACT/RES DGA 2240 CS GAMMA = RES DGA 2250 DGAMMA = RES DGA 2260 120 RETURN DGA 2270 C ---------- LAST CARD OF GAMMA ---------- DGA 2280 END DGA 2290 subroutine r8vec_print ( n, a, title ) c*********************************************************************72 c cc R8VEC_PRINT prints an R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of components of the vector. c c Input, double precision A(N), the vector to be printed. c c Input, character * ( * ) TITLE, a title. c implicit none integer n double precision a(n) integer i integer s_len_trim character ( len = * ) title integer title_length write ( *, '(a)' ) ' ' write ( *, '(a)' ) title write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,g16.8)' ) i, a(i) end do return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end