# To unbundle, sh this file echo iqd.f 1>&2 cat >iqd.f <<'End of iqd.f' 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 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 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 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 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 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 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 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 C SUBROUTINE CLIQFS(NT,T,WTS,KIND,ALPHA,BETA, 1LO,NWF,WF,NIWF,IWF,IER) 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 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 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 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 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 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 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 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 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 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 DO20J=1,NSTAR WF(J)=A(J) IF(NM.LT.1) GOTO 20 DO10I=1,NM WF(J)=WF(J)*(PHI(J)-XK(I)) 10 CONTINUE 20 CONTINUE C COMPUTE Y-HAT DO40I=1,M SUM=ZERO DO30J=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 DO70I=1,NM TMP=V-XK(I) IF(L.EQ.1) GOTO 60 LAST=MIN(L,I+1) DO50JR=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 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 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 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 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 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 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 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 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 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 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 (ESSENTIALY) 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 110 RETURN END SUBROUTINE MACHEP (X) 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---------------------------------------------------------------------- DGA 30 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 End of iqd.f echo iqtest.d 1>&2 cat >iqtest.d <<'End of iqtest.d' ---------------------------------------- TESTING CIQFS 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 1 (-1,1) ONE LEGENDRE MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 0.95105651629515357D+00 2 0.22240110861588534D+00 -0.73134471884532099D-02 2 0.58778525229247312D+00 2 0.48363063741586058D+00 -0.17871860197559908D-01 3 -0.22034386889519082D-16 2 0.58793650793650781D+00 0.86736173798840349D-17 4 -0.58778525229247315D+00 2 0.48363063741586064D+00 0.17871860197559929D-01 5 -0.95105651629515357D+00 2 0.22240110861588531D+00 0.73134471884532081D-02 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.104D-16 0.104D-16 MAXIMUM : 0.305D-15 0.102D-15 WEIGHTS RATIO 0.667D+00 ERROR FOR 10-TH POWER 0.387D-02 ERROR CONSTANT 0.107D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2000000000D+01 0.2000000000D+01 0.305D-15 0.102D-15 2 0.0000000000D+00 -0.1387778781D-16 0.139D-16 0.139D-16 3 0.6666666667D+00 0.6666666667D+00 0.278D-16 0.167D-16 4 0.0000000000D+00 0.2081668171D-16 -0.208D-16 -0.208D-16 5 0.4000000000D+00 0.4000000000D+00 -0.347D-16 -0.248D-16 6 0.0000000000D+00 0.1734723476D-16 -0.173D-16 -0.173D-16 7 0.2857142857D+00 0.2857142857D+00 -0.555D-16 -0.432D-16 8 0.0000000000D+00 0.1387778781D-16 -0.139D-16 -0.139D-16 9 0.2222222222D+00 0.2222222222D+00 -0.590D-16 -0.483D-16 10 0.0000000000D+00 0.1040834086D-16 -0.104D-16 -0.104D-16 11 0.1818181818D+00 0.1779513889D+00 0.387D-02 0.327D-02 12 0.0000000000D+00 0.3469446952D-17 -0.347D-17 -0.347D-17 13 0.1538461538D+00 0.1429191468D+00 0.109D-01 0.947D-02 RETURN FROM CIQFS: ERROR INDICATOR IER= 0 ---------------------------------------- TESTING CIQF,CIQFS,CGQF AND CGQFS WITH ALL CLASSICAL WEIGHT FUNCTIONS KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 1 (A,B) ONE LEGENDRE PARAMETERS A -0.50000 B 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.38272480742333004D+00 1 0.29615860632023627D+00 2 0.76913362367896129D-01 1 0.59828583812420807D+00 3 0.74999999999999999D+00 1 0.71111111111111097D+00 4 0.14230866376321039D+01 1 0.59828583812420795D+00 5 0.18827248074233300D+01 1 0.29615860632023639D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.416D-16 0.241D-16 MAXIMUM : 0.389D-15 0.375D-15 WEIGHTS RATIO 0.714D+00 ERROR FOR 10-TH POWER 0.341D-01 ERROR CONSTANT 0.941D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2500000000D+01 0.2500000000D+01 0.333D-15 0.952D-16 2 0.0000000000D+00 0.4163336342D-16 -0.416D-16 -0.416D-16 3 0.1302083333D+01 0.1302083333D+01 0.555D-16 0.241D-16 4 0.0000000000D+00 0.1457167720D-15 -0.146D-15 -0.146D-15 5 0.1220703125D+01 0.1220703125D+01 -0.555D-16 -0.250D-16 6 0.0000000000D+00 0.2081668171D-15 -0.208D-15 -0.208D-15 7 0.1362391881D+01 0.1362391881D+01 -0.222D-15 -0.940D-16 8 0.0000000000D+00 0.2636779683D-15 -0.264D-15 -0.264D-15 9 0.1655684577D+01 0.1655684577D+01 -0.389D-15 -0.146D-15 10 0.0000000000D+00 0.3747002708D-15 -0.375D-15 -0.375D-15 11 0.2116642215D+01 0.2082511426D+01 0.341D-01 0.110D-01 12 0.0000000000D+00 0.4440892099D-15 -0.444D-15 -0.444D-15 13 0.2798445236D+01 0.2653042970D+01 0.145D+00 0.383D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 1 (A,B) ONE LEGENDRE PARAMETERS A -0.50000 B 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.38272480742333004D+00 2 0.29615860632023633D+00 0.51375328699628709D-17 2 0.76913362367896129D-01 2 0.59828583812420802D+00 0.10378604887385233D-16 3 0.74999999999999999D+00 2 0.71111111111111097D+00 0.70271193197689065D-17 4 0.14230866376321039D+01 2 0.59828583812420796D+00 -0.10378604887385234D-16 5 0.18827248074233300D+01 2 0.29615860632023636D+00 -0.10517194876048111D-33 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.694D-17 0.694D-17 MAXIMUM : 0.361D-15 0.222D-15 WEIGHTS RATIO 0.714D+00 ERROR FOR 10-TH POWER 0.341D-01 ERROR CONSTANT 0.941D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2500000000D+01 0.2500000000D+01 0.333D-15 0.952D-16 2 0.0000000000D+00 0.6938893904D-17 -0.694D-17 -0.694D-17 3 0.1302083333D+01 0.1302083333D+01 0.833D-16 0.362D-16 4 0.0000000000D+00 0.6938893904D-16 -0.694D-16 -0.694D-16 5 0.1220703125D+01 0.1220703125D+01 -0.278D-16 -0.125D-16 6 0.0000000000D+00 0.1110223025D-15 -0.111D-15 -0.111D-15 7 0.1362391881D+01 0.1362391881D+01 -0.194D-15 -0.822D-16 8 0.0000000000D+00 0.1526556659D-15 -0.153D-15 -0.153D-15 9 0.1655684577D+01 0.1655684577D+01 -0.361D-15 -0.136D-15 10 0.0000000000D+00 0.2220446049D-15 -0.222D-15 -0.222D-15 11 0.2116642215D+01 0.2082511426D+01 0.341D-01 0.110D-01 12 0.0000000000D+00 0.3330669074D-15 -0.333D-15 -0.333D-15 13 0.2798445236D+01 0.2653042970D+01 0.145D+00 0.383D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 2 (A,B) ((B-X)*(X-A))**(-HALF) CHEBYSHEV PARAMETERS A -0.50000 B 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.43882064536894194D+00 1 0.62831853071795850D+00 2 0.15268434634408592D-01 1 0.62831853071795866D+00 3 0.75000000000000000D+00 1 0.62831853071795835D+00 4 0.14847315653655914D+01 1 0.62831853071795842D+00 5 0.19388206453689420D+01 1 0.62831853071795848D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.278D-16 0.278D-16 MAXIMUM : 0.167D-14 0.444D-15 WEIGHTS RATIO 0.759D+00 ERROR FOR 10-TH POWER 0.571D-01 ERROR CONSTANT 0.157D-07 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3141592654D+01 0.3141592654D+01 0.777D-15 0.188D-15 2 0.0000000000D+00 -0.1804112415D-15 0.180D-15 0.180D-15 3 0.2454369261D+01 0.2454369261D+01 0.555D-15 0.161D-15 4 0.0000000000D+00 -0.2775557562D-16 0.278D-16 0.278D-16 5 0.2876213977D+01 0.2876213977D+01 0.722D-15 0.186D-15 6 0.0000000000D+00 0.1110223025D-15 -0.111D-15 -0.111D-15 7 0.3745070283D+01 0.3745070283D+01 0.111D-14 0.234D-15 8 0.0000000000D+00 0.2775557562D-15 -0.278D-15 -0.278D-15 9 0.5120213277D+01 0.5120213277D+01 0.167D-14 0.272D-15 10 0.0000000000D+00 0.4440892099D-15 -0.444D-15 -0.444D-15 11 0.7200299921D+01 0.7143154684D+01 0.571D-01 0.697D-02 12 0.0000000000D+00 0.7771561172D-15 -0.777D-15 -0.777D-15 13 0.1031292957D+02 0.1004506127D+02 0.268D+00 0.237D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 2 (A,B) ((B-X)*(X-A))**(-HALF) CHEBYSHEV PARAMETERS A -0.50000 B 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.43882064536894194D+00 2 0.62831853071795850D+00 0.15624107072989836D-33 2 0.15268434634408592D-01 2 0.62831853071795868D+00 0.10899589056276976D-16 3 0.75000000000000000D+00 2 0.62831853071795843D+00 -0.35666697511613738D-17 4 0.14847315653655914D+01 2 0.62831853071795848D+00 -0.10899589056276969D-16 5 0.19388206453689420D+01 2 0.62831853071795851D+00 -0.10899589056276971D-16 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.000D+00 0.000D+00 MAXIMUM : 0.189D-14 0.308D-15 WEIGHTS RATIO 0.759D+00 ERROR FOR 10-TH POWER 0.571D-01 ERROR CONSTANT 0.157D-07 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3141592654D+01 0.3141592654D+01 0.611D-15 0.147D-15 2 0.0000000000D+00 -0.1110223025D-15 0.111D-15 0.111D-15 3 0.2454369261D+01 0.2454369261D+01 0.555D-15 0.161D-15 4 0.0000000000D+00 0.0000000000D+00 0.000D+00 0.000D+00 5 0.2876213977D+01 0.2876213977D+01 0.777D-15 0.200D-15 6 0.0000000000D+00 0.5551115123D-16 -0.555D-16 -0.555D-16 7 0.3745070283D+01 0.3745070283D+01 0.122D-14 0.257D-15 8 0.0000000000D+00 0.1110223025D-15 -0.111D-15 -0.111D-15 9 0.5120213277D+01 0.5120213277D+01 0.189D-14 0.308D-15 10 0.0000000000D+00 0.2220446049D-15 -0.222D-15 -0.222D-15 11 0.7200299921D+01 0.7143154684D+01 0.571D-01 0.697D-02 12 0.0000000000D+00 0.3330669074D-15 -0.333D-15 -0.333D-15 13 0.1031292957D+02 0.1004506127D+02 0.268D+00 0.237D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 3 (A,B) ((B-X)*(X-A))**ALPHA GEGENBAUER PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.33253175473054836D+00 1 0.20453077171808547D+00 2 0.12500000000000004D+00 1 0.61359231515425650D+00 3 0.75000000000000000D+00 1 0.81812308687234188D+00 4 0.13750000000000000D+01 1 0.61359231515425636D+00 5 0.18325317547305483D+01 1 0.20453077171808545D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.139D-16 0.771D-17 MAXIMUM : 0.333D-15 0.146D-15 WEIGHTS RATIO 0.711D+00 ERROR FOR 10-TH POWER 0.223D-01 ERROR CONSTANT 0.615D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2454369261D+01 0.2454369261D+01 0.333D-15 0.964D-16 2 0.0000000000D+00 -0.8673617380D-16 0.867D-16 0.867D-16 3 0.9587379924D+00 0.9587379924D+00 0.125D-15 0.638D-16 4 0.0000000000D+00 -0.5551115123D-16 0.555D-16 0.555D-16 5 0.7490140566D+00 0.7490140566D+00 0.833D-16 0.476D-16 6 0.0000000000D+00 -0.6245004514D-16 0.625D-16 0.625D-16 7 0.7314590396D+00 0.7314590396D+00 0.139D-16 0.802D-17 8 0.0000000000D+00 -0.1040834086D-15 0.104D-15 0.104D-15 9 0.8000333246D+00 0.8000333246D+00 -0.139D-16 -0.771D-17 10 0.0000000000D+00 -0.1457167720D-15 0.146D-15 0.146D-15 11 0.9375390523D+00 0.9152166939D+00 0.223D-01 0.115D-01 12 0.0000000000D+00 -0.2012279232D-15 0.201D-15 0.201D-15 13 0.1150996604D+01 0.1063799892D+01 0.872D-01 0.405D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 3 (A,B) ((B-X)*(X-A))**ALPHA GEGENBAUER PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.33253175473054836D+00 2 0.20453077171808549D+00 0.51690110349341779D-34 2 0.12500000000000004D+00 2 0.61359231515425650D+00 0.53220649688852420D-17 3 0.75000000000000000D+00 2 0.81812308687234199D+00 0.71948497126933839D-18 4 0.13750000000000000D+01 2 0.61359231515425630D+00 -0.10644129937770478D-16 5 0.18325317547305483D+01 2 0.20453077171808546D+00 0.35480433125901607D-17 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.416D-16 0.238D-16 MAXIMUM : 0.222D-15 0.135D-15 WEIGHTS RATIO 0.711D+00 ERROR FOR 10-TH POWER 0.223D-01 ERROR CONSTANT 0.615D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2454369261D+01 0.2454369261D+01 0.222D-15 0.643D-16 2 0.0000000000D+00 -0.1353084311D-15 0.135D-15 0.135D-15 3 0.9587379924D+00 0.9587379924D+00 0.125D-15 0.638D-16 4 0.0000000000D+00 -0.6938893904D-16 0.694D-16 0.694D-16 5 0.7490140566D+00 0.7490140566D+00 0.416D-16 0.238D-16 6 0.0000000000D+00 -0.5551115123D-16 0.555D-16 0.555D-16 7 0.7314590396D+00 0.7314590396D+00 -0.416D-16 -0.240D-16 8 0.0000000000D+00 -0.8326672685D-16 0.833D-16 0.833D-16 9 0.8000333246D+00 0.8000333246D+00 -0.125D-15 -0.694D-16 10 0.0000000000D+00 -0.9714451465D-16 0.971D-16 0.971D-16 11 0.9375390523D+00 0.9152166939D+00 0.223D-01 0.115D-01 12 0.0000000000D+00 -0.1318389842D-15 0.132D-15 0.132D-15 13 0.1150996604D+01 0.1063799892D+01 0.872D-01 0.405D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 4 (A,B) (B-X)**ALPHA*(X-A)**BETA JACOBI PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 BETA 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.15303633066793587D+00 1 0.76617129890421453D-01 2 0.35753707671627774D+00 1 0.53659067985975822D+00 3 0.94535638653559417D+00 1 0.12551248538055703D+01 4 0.14858880884829689D+01 1 0.13472898724319172D+01 5 0.18642547789330951D+01 1 0.54899372611754609D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.000D+00 0.000D+00 MAXIMUM : 0.333D-15 0.775D-16 WEIGHTS RATIO 0.790D+00 ERROR FOR 10-TH POWER 0.146D-01 ERROR CONSTANT 0.402D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3764616262D+01 0.3764616262D+01 0.333D-15 0.699D-16 2 0.1568590109D+01 0.1568590109D+01 0.000D+00 0.000D+00 3 0.1604239884D+01 0.1604239884D+01 0.278D-16 0.107D-16 4 0.1216891365D+01 0.1216891365D+01 0.278D-16 0.125D-16 5 0.1306872769D+01 0.1306872769D+01 0.555D-16 0.241D-16 6 0.1183053821D+01 0.1183053821D+01 0.278D-16 0.127D-16 7 0.1308228360D+01 0.1308228360D+01 0.555D-16 0.240D-16 8 0.1289910261D+01 0.1289910261D+01 0.111D-15 0.485D-16 9 0.1454550385D+01 0.1454550385D+01 0.111D-15 0.452D-16 10 0.1508092819D+01 0.1508092819D+01 0.194D-15 0.775D-16 11 0.1724613987D+01 0.1710038800D+01 0.146D-01 0.535D-02 12 0.1848110450D+01 0.1825870725D+01 0.222D-01 0.781D-02 13 0.2135936129D+01 0.2067375141D+01 0.686D-01 0.219D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 4 (A,B) (B-X)**ALPHA*(X-A)**BETA JACOBI PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 BETA 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.15303633066793587D+00 2 0.76617129890421451D-01 0.32590860751743046D-34 2 0.35753707671627774D+00 2 0.53659067985975821D+00 0.46541822467153897D-17 3 0.94535638653559417D+00 2 0.12551248538055703D+01 0.17784281776530149D-32 4 0.14858880884829689D+01 2 0.13472898724319171D+01 0.33436908629561223D-32 5 0.18642547789330951D+01 2 0.54899372611754609D+00 -0.95235230486008862D-17 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.278D-16 0.108D-16 MAXIMUM : 0.389D-15 0.155D-15 WEIGHTS RATIO 0.790D+00 ERROR FOR 10-TH POWER 0.146D-01 ERROR CONSTANT 0.402D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3764616262D+01 0.3764616262D+01 0.389D-15 0.816D-16 2 0.1568590109D+01 0.1568590109D+01 0.278D-16 0.108D-16 3 0.1604239884D+01 0.1604239884D+01 0.833D-16 0.320D-16 4 0.1216891365D+01 0.1216891365D+01 0.833D-16 0.376D-16 5 0.1306872769D+01 0.1306872769D+01 0.111D-15 0.481D-16 6 0.1183053821D+01 0.1183053821D+01 0.111D-15 0.509D-16 7 0.1308228360D+01 0.1308228360D+01 0.139D-15 0.601D-16 8 0.1289910261D+01 0.1289910261D+01 0.250D-15 0.109D-15 9 0.1454550385D+01 0.1454550385D+01 0.278D-15 0.113D-15 10 0.1508092819D+01 0.1508092819D+01 0.389D-15 0.155D-15 11 0.1724613987D+01 0.1710038800D+01 0.146D-01 0.535D-02 12 0.1848110450D+01 0.1825870725D+01 0.222D-01 0.781D-02 13 0.2135936129D+01 0.2067375141D+01 0.686D-01 0.219D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 5 (A,INF) (X-A)**ALPHA*EXP(-B*(X-A)) GEN LAGUERRE PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.28430059642607425D+00 1 0.13097405507334787D+00 2 0.37987684921184824D+00 1 0.14587060425199263D+00 3 0.15522326814141575D+01 1 0.34570386911402214D-01 4 0.33733518897712785D+01 1 0.18997892129372685D-02 5 0.62288391760287902D+01 1 0.13698879195062384D-04 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.104D-16 0.843D-17 MAXIMUM : 0.995D-13 0.156D-15 WEIGHTS RATIO 0.239D+00 ERROR FOR 10-TH POWER 0.119D+02 ERROR CONSTANT 0.329D-05 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3133285343D+00 0.3133285343D+00 0.139D-16 0.106D-16 2 0.2349964007D+00 0.2349964007D+00 0.104D-16 0.843D-17 3 0.2937455009D+00 0.2937455009D+00 0.208D-16 0.161D-16 4 0.5140546266D+00 0.5140546266D+00 0.555D-16 0.367D-16 5 0.1156622910D+01 0.1156622910D+01 0.167D-15 0.772D-16 6 0.3180713002D+01 0.3180713002D+01 0.500D-15 0.120D-15 7 0.1033731726D+02 0.1033731726D+02 0.155D-14 0.137D-15 8 0.3876493972D+02 0.3876493972D+02 0.622D-14 0.156D-15 9 0.1647509938D+03 0.1647509938D+03 0.249D-13 0.150D-15 10 0.7825672205D+03 0.7825672205D+03 0.995D-13 0.127D-15 11 0.4108477908D+04 0.4096550234D+04 0.119D+02 0.290D-02 12 0.2362374797D+05 0.2322715282D+05 0.397D+03 0.168D-01 13 0.1476484248D+06 0.1398805273D+06 0.777D+04 0.526D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 5 (A,INF) (X-A)**ALPHA*EXP(-B*(X-A)) GEN LAGUERRE PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.28430059642607425D+00 2 0.13097405507334787D+00 0.45440753615923155D-18 2 0.37987684921184824D+00 2 0.14587060425199263D+00 -0.10783931190249757D-34 3 0.15522326814141575D+01 2 0.34570386911402215D-01 -0.21557758188527963D-35 4 0.33733518897712785D+01 2 0.18997892129372686D-02 -0.11226942266757488D-36 5 0.62288391760287902D+01 2 0.13698879195062384D-04 0.00000000000000000D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.173D-16 0.140D-16 MAXIMUM : 0.853D-13 0.118D-15 WEIGHTS RATIO 0.239D+00 ERROR FOR 10-TH POWER 0.119D+02 ERROR CONSTANT 0.329D-05 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.3133285343D+00 0.3133285343D+00 0.208D-16 0.159D-16 2 0.2349964007D+00 0.2349964007D+00 0.173D-16 0.140D-16 3 0.2937455009D+00 0.2937455009D+00 0.208D-16 0.161D-16 4 0.5140546266D+00 0.5140546266D+00 0.555D-16 0.367D-16 5 0.1156622910D+01 0.1156622910D+01 0.139D-15 0.643D-16 6 0.3180713002D+01 0.3180713002D+01 0.389D-15 0.929D-16 7 0.1033731726D+02 0.1033731726D+02 0.133D-14 0.118D-15 8 0.3876493972D+02 0.3876493972D+02 0.444D-14 0.112D-15 9 0.1647509938D+03 0.1647509938D+03 0.178D-13 0.107D-15 10 0.7825672205D+03 0.7825672205D+03 0.853D-13 0.109D-15 11 0.4108477908D+04 0.4096550234D+04 0.119D+02 0.290D-02 12 0.2362374797D+05 0.2322715282D+05 0.397D+03 0.168D-01 13 0.1476484248D+06 0.1398805273D+06 0.777D+04 0.526D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 6 (-INF,INF) ABS(X-A)**ALFA*EXP(-B*(X-A)**2) GEN HERMITE PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.19846400902538129D+01 1 0.12302854647083805D-01 2 -0.12388124270822399D+01 1 0.20061059263754047D+00 3 -0.49999999999999999D+00 1 0.30281023613813227D+00 4 0.23881242708223986D+00 1 0.20061059263754044D+00 5 0.98464009025381294D+00 1 0.12302854647083807D-01 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.347D-17 0.347D-17 MAXIMUM : 0.236D-15 0.146D-15 WEIGHTS RATIO 0.422D+00 ERROR FOR 10-TH POWER 0.164D+00 ERROR CONSTANT 0.453D-07 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.7286371307D+00 0.7286371307D+00 0.694D-16 0.401D-16 2 0.0000000000D+00 -0.1431146868D-16 0.143D-16 0.143D-16 3 0.2732389240D+00 0.2732389240D+00 0.347D-16 0.272D-16 4 0.0000000000D+00 -0.1214306433D-16 0.121D-16 0.121D-16 5 0.2390840585D+00 0.2390840585D+00 0.520D-16 0.420D-16 6 0.0000000000D+00 -0.3469446952D-17 0.347D-17 0.347D-17 7 0.3287405805D+00 0.3287405805D+00 0.104D-15 0.783D-16 8 0.0000000000D+00 0.1387778781D-16 -0.139D-16 -0.139D-16 9 0.6163885884D+00 0.6163885884D+00 0.236D-15 0.146D-15 10 0.0000000000D+00 0.4857225733D-16 -0.486D-16 -0.486D-16 11 0.1463922897D+01 0.1299552607D+01 0.164D+00 0.667D-01 12 0.0000000000D+00 0.9714451465D-16 -0.971D-16 -0.971D-16 13 0.4208778330D+01 0.2832177149D+01 0.138D+01 0.264D+00 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 6 (-INF,INF) ABS(X-A)**ALFA*EXP(-B*(X-A)**2) GEN HERMITE PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.19846400902538129D+01 2 0.12302854647083803D-01 0.48291548254172280D-18 2 -0.12388124270822399D+01 2 0.20061059263754046D+00 0.39372147329035634D-17 3 -0.49999999999999999D+00 2 0.30281023613813223D+00 0.63813859671973973D-18 4 0.23881242708223986D+00 2 0.20061059263754042D+00 -0.19211226557309852D-33 5 0.98464009025381294D+00 2 0.12302854647083805D-01 -0.85397410072820073D-35 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.260D-17 0.260D-17 MAXIMUM : 0.389D-15 0.240D-15 WEIGHTS RATIO 0.422D+00 ERROR FOR 10-TH POWER 0.164D+00 ERROR CONSTANT 0.453D-07 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.7286371307D+00 0.7286371307D+00 0.125D-15 0.723D-16 2 0.0000000000D+00 -0.1951563910D-16 0.195D-16 0.195D-16 3 0.2732389240D+00 0.2732389240D+00 0.694D-16 0.545D-16 4 0.0000000000D+00 -0.2602085214D-17 0.260D-17 0.260D-17 5 0.2390840585D+00 0.2390840585D+00 0.937D-16 0.756D-16 6 0.0000000000D+00 0.2255140519D-16 -0.226D-16 -0.226D-16 7 0.3287405805D+00 0.3287405805D+00 0.173D-15 0.131D-15 8 0.0000000000D+00 0.6938893904D-16 -0.694D-16 -0.694D-16 9 0.6163885884D+00 0.6163885884D+00 0.389D-15 0.240D-15 10 0.0000000000D+00 0.1873501354D-15 -0.187D-15 -0.187D-15 11 0.1463922897D+01 0.1299552607D+01 0.164D+00 0.667D-01 12 0.0000000000D+00 0.4857225733D-15 -0.486D-15 -0.486D-15 13 0.4208778330D+01 0.2832177149D+01 0.138D+01 0.264D+00 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 7 (A,B) ABS(X-(A+B)/TWO)**ALFA EXPONENTIAL PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.39148162917747822D+00 1 0.29332908318369346D+00 2 0.38501428978160557D-01 1 0.47745248689814077D+00 3 0.75000000000000006D+00 1 0.32182684108615619D+00 4 0.14614985710218395D+01 1 0.47745248689814066D+00 5 0.18914816291774782D+01 1 0.29332908318369348D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.416D-16 0.416D-16 MAXIMUM : 0.103D-14 0.375D-15 WEIGHTS RATIO 0.651D+00 ERROR FOR 10-TH POWER 0.285D-01 ERROR CONSTANT 0.786D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.1863389981D+01 0.1863389981D+01 0.194D-15 0.679D-16 2 0.0000000000D+00 -0.4163336342D-16 0.416D-16 0.416D-16 3 0.1247805791D+01 0.1247805791D+01 0.222D-15 0.988D-16 4 0.0000000000D+00 -0.5551115123D-16 0.555D-16 0.555D-16 5 0.1240715985D+01 0.1240715985D+01 0.416D-15 0.186D-15 6 0.0000000000D+00 -0.9714451465D-16 0.971D-16 0.971D-16 7 0.1421653733D+01 0.1421653733D+01 0.694D-15 0.287D-15 8 0.0000000000D+00 -0.1942890293D-15 0.194D-15 0.194D-15 9 0.1753684704D+01 0.1753684704D+01 0.103D-14 0.373D-15 10 0.0000000000D+00 -0.3747002708D-15 0.375D-15 0.375D-15 11 0.2263587593D+01 0.2235050644D+01 0.285D-01 0.874D-02 12 0.0000000000D+00 -0.5828670879D-15 0.583D-15 0.583D-15 13 0.3012877005D+01 0.2886932684D+01 0.126D+00 0.314D-01 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 7 (A,B) ABS(X-(A+B)/TWO)**ALFA EXPONENTIAL PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.39148162917747822D+00 2 0.29332908318369348D+00 0.50529291023025702D-34 2 0.38501428978160557D-01 2 0.47745248689814078D+00 0.82824803768571371D-17 3 0.75000000000000006D+00 2 0.32182684108615620D+00 0.57575818627259635D-18 4 0.14614985710218395D+01 2 0.47745248689814064D+00 -0.10193662741455050D-33 5 0.18914816291774782D+01 2 0.29332908318369348D+00 -0.11795232945122243D-34 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.902D-16 0.582D-16 MAXIMUM : 0.971D-15 0.444D-15 WEIGHTS RATIO 0.651D+00 ERROR FOR 10-TH POWER 0.285D-01 ERROR CONSTANT 0.786D-08 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.1863389981D+01 0.1863389981D+01 0.167D-15 0.582D-16 2 0.0000000000D+00 -0.9020562075D-16 0.902D-16 0.902D-16 3 0.1247805791D+01 0.1247805791D+01 0.222D-15 0.988D-16 4 0.0000000000D+00 -0.9020562075D-16 0.902D-16 0.902D-16 5 0.1240715985D+01 0.1240715985D+01 0.389D-15 0.173D-15 6 0.0000000000D+00 -0.1387778781D-15 0.139D-15 0.139D-15 7 0.1421653733D+01 0.1421653733D+01 0.638D-15 0.264D-15 8 0.0000000000D+00 -0.2498001805D-15 0.250D-15 0.250D-15 9 0.1753684704D+01 0.1753684704D+01 0.971D-15 0.353D-15 10 0.0000000000D+00 -0.4440892099D-15 0.444D-15 0.444D-15 11 0.2263587593D+01 0.2235050644D+01 0.285D-01 0.874D-02 12 0.0000000000D+00 -0.6661338148D-15 0.666D-15 0.666D-15 13 0.3012877005D+01 0.2886932684D+01 0.126D+00 0.314D-01 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 KNOTS AND WEIGHTS OF GAUSS QF COMPUTED BY CGQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 8 (A,INF) (X-A)**ALFA*(B+X)**BETA RATIONAL PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 BETA -16.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.43541078852037897D+00 1 0.26216877717069183D-04 2 -0.21664377601574750D+00 1 0.16250543349707103D-04 3 0.25596297684363936D+00 1 0.12927396984510148D-05 4 0.12864478502358691D+01 1 0.11152074844159659D-07 5 0.41096437374566179D+01 1 0.28799054822892455D-11 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.265D-22 0.265D-22 MAXIMUM : 0.254D-20 0.254D-20 WEIGHTS RATIO 0.438D-04 ERROR FOR 10-TH POWER 0.825D-06 ERROR CONSTANT 0.227D-12 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.4377131572D-04 0.4377131572D-04 0.254D-20 0.254D-20 2 0.7295219287D-05 0.7295219287D-05 0.318D-21 0.318D-21 3 0.2188565786D-05 0.2188565786D-05 0.529D-22 0.529D-22 4 0.9991278588D-06 0.9991278588D-06 0.265D-22 0.265D-22 5 0.6422964807D-06 0.6422964807D-06 0.397D-22 0.397D-22 6 0.5577837858D-06 0.5577837858D-06 0.662D-22 0.662D-22 7 0.6398108132D-06 0.6398108132D-06 0.159D-21 0.159D-21 8 0.9597162198D-06 0.9597162198D-06 0.371D-21 0.371D-21 9 0.1882520277D-05 0.1882520277D-05 0.900D-21 0.900D-21 10 0.4877438900D-05 0.4877438900D-05 0.244D-20 0.244D-20 11 0.1707103615D-04 0.1624615493D-04 0.825D-06 0.825D-06 12 0.8413582103D-04 0.6416191151D-04 0.200D-04 0.200D-04 13 0.6310186577D-03 0.2769134667D-03 0.354D-03 0.354D-03 RETURN FROM CGQF(S): ERROR INDICATOR IER= 0 WEIGHTS OF GAUSS QF COMPUTED FROM THE KNOTS BY CIQF(S) 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 8 (A,INF) (X-A)**ALFA*(B+X)**BETA RATIONAL PARAMETERS A -0.50000 B 2.00000 ALPHA 0.50000 BETA -16.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 -0.43541078852037897D+00 2 0.26216877717069183D-04 0.68218549863919726D-22 2 -0.21664377601574750D+00 2 0.16250543349707102D-04 0.35129201900358245D-38 3 0.25596297684363936D+00 2 0.12927396984510146D-05 -0.26910550838761838D-22 4 0.12864478502358691D+01 2 0.11152074844159657D-07 -0.46429838491233957D-24 5 0.41096437374566179D+01 2 0.28799054822892453D-11 0.00000000000000000D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 10 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.185D-21 0.185D-21 MAXIMUM : 0.424D-20 0.423D-20 WEIGHTS RATIO 0.438D-04 ERROR FOR 10-TH POWER 0.825D-06 ERROR CONSTANT 0.227D-12 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.4377131572D-04 0.4377131572D-04 0.424D-20 0.423D-20 2 0.7295219287D-05 0.7295219287D-05 0.741D-21 0.741D-21 3 0.2188565786D-05 0.2188565786D-05 0.318D-21 0.318D-21 4 0.9991278588D-06 0.9991278588D-06 0.212D-21 0.212D-21 5 0.6422964807D-06 0.6422964807D-06 0.185D-21 0.185D-21 6 0.5577837858D-06 0.5577837858D-06 0.212D-21 0.212D-21 7 0.6398108132D-06 0.6398108132D-06 0.344D-21 0.344D-21 8 0.9597162198D-06 0.9597162198D-06 0.635D-21 0.635D-21 9 0.1882520277D-05 0.1882520277D-05 0.132D-20 0.132D-20 10 0.4877438900D-05 0.4877438900D-05 0.339D-20 0.339D-20 11 0.1707103615D-04 0.1624615493D-04 0.825D-06 0.825D-06 12 0.8413582103D-04 0.6416191151D-04 0.200D-04 0.200D-04 13 0.6310186577D-03 0.2769134667D-03 0.354D-03 0.354D-03 RETURN FROM CIQF(S): ERROR INDICATOR IER= 0 ---------------------------------------- TESTING CEIQFS RETURN FROM CEIQFS: ERROR INDICATOR IER= 0 INTERGRAL OF SIN(X) ON [-1,1] BY FEJER TYPE RULE WITH 5 POINTS OF MULTIPLICITY 2 QUADRATURE FORMULA: -1.7564075194265172E-17 EXACT VALUE : 0.0000000000000000E+00 ---------------------------------------- TESTING CEIQF RETURN FROM CEIQF: ERROR INDICATOR IER= 0 INTERGRAL OF SIN(X) FROM -0.5000000000000000 TO 2.000000000000000 BY FEJER TYPE RULE WITH 5 POINTS OF MULTIPLICITY 2 QUADRATURE FORMULA: 1.293729406614736 EXACT VALUE : 1.293729398437515 ---------------------------------------- TESTING CLIQFS 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 1 (-1,1) ONE LEGENDRE MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 0.95105651629515357D+00 1 0.16778122846668345D+00 2 0.58778525229247312D+00 1 0.52555210486664980D+00 3 -0.22034386889519082D-16 1 0.61333333333333317D+00 4 -0.58778525229247315D+00 1 0.52555210486664973D+00 5 -0.95105651629515357D+00 1 0.16778122846668345D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 5 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.347D-17 0.347D-17 MAXIMUM : 0.361D-15 0.120D-15 WEIGHTS RATIO 0.667D+00 ERROR FOR 5-TH POWER -0.347D-17 ERROR CONSTANT -0.289D-19 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2000000000D+01 0.2000000000D+01 0.361D-15 0.120D-15 2 0.0000000000D+00 0.1734723476D-16 -0.173D-16 -0.173D-16 3 0.6666666667D+00 0.6666666667D+00 0.111D-15 0.666D-16 4 0.0000000000D+00 0.3469446952D-17 -0.347D-17 -0.347D-17 5 0.4000000000D+00 0.4000000000D+00 0.763D-16 0.545D-16 6 0.0000000000D+00 0.3469446952D-17 -0.347D-17 -0.347D-17 7 0.2857142857D+00 0.2916666667D+00 -0.595D-02 -0.463D-02 8 0.0000000000D+00 0.0000000000D+00 0.000D+00 0.000D+00 RETURN FROM CLIQFS. IER= 0 ---------------------------------------- TESTING CLIQF AND EIQFS 1 INTERPOLATORY QUADRATURE FORMULA TYPE INTERVAL WEIGHT FUNCTION NAME 1 (A,B) ONE LEGENDRE PARAMETERS A -0.50000 B 2.00000 MACHINE PRECISION 0.1D-16 KNOTS MULT WEIGHTS 1 0.19388206453689420D+01 1 0.20972653558335429D+00 2 0.14847315653655914D+01 1 0.65694013108331227D+00 3 0.74999999999999997D+00 1 0.76666666666666645D+00 4 0.15268434634408565D-01 1 0.65694013108331215D+00 5 -0.43882064536894197D+00 1 0.20972653558335429D+00 COMPARISON OF MOMENTS ORDER OF PRECISION 5 ERRORS : ABSOLUTE RELATIVE ---------+------------------------- MINIMUM : 0.486D-16 0.486D-16 MAXIMUM : 0.500D-15 0.143D-15 WEIGHTS RATIO 0.714D+00 ERROR FOR 5-TH POWER -0.694D-17 ERROR CONSTANT -0.578D-19 MOMENTS: TRUE FROM Q.F. ERROR RELATIVE 1 0.2500000000D+01 0.2500000000D+01 0.500D-15 0.143D-15 2 0.0000000000D+00 0.5204170428D-16 -0.520D-16 -0.520D-16 3 0.1302083333D+01 0.1302083333D+01 0.278D-15 0.121D-15 4 0.0000000000D+00 0.4857225733D-16 -0.486D-16 -0.486D-16 5 0.1220703125D+01 0.1220703125D+01 0.278D-15 0.125D-15 6 0.0000000000D+00 0.6938893904D-17 -0.694D-17 -0.694D-17 7 0.1362391881D+01 0.1390775045D+01 -0.284D-01 -0.120D-01 8 0.0000000000D+00 -0.1387778781D-16 0.139D-16 0.139D-16 RETURN FROM CLIQF: ERROR INDICATOR IER= 0 RETURN FROM EIQFS: ERROR INDICATOR IER= 0 INTERGRAL OF SIN(X) FROM -0.5000000000000000 TO 2.000000000000000 BY FEJER TYPE RULE WITH 5 POINTS OF MULTIPLICITY ONE QUADRATURE FORMULA: 1.293704657106341 EXACT VALUE : 1.293729398437515 ---------------------------------------- TESTING CEGQF RETURN FROM CEGQF: ERROR INDICATOR IER= 0 INTERGRAL OF X*SIN(X) FROM 1.000000000000000 TO 4.000000000000000 BY GAUSS-EXPONENTIAL RULE WITH 12 POINTS QUADRATURE FORMULA: 0.6786430913896829 EXACT VALUE : 0.6786430913896827 TESTING CEGQFS RETURN FROM CEGQFS: ERROR INDICATOR IER= 0 INTERGRAL OF X*SIN(X) FROM -1.000000000000000 TO 1.000000000000000 BY GAUSS-EXPONENTIAL RULE WITH 12 POINTS QUADRATURE FORMULA: -1.7694179454963432E-16 EXACT VALUE : 0.0000000000000000E+00 End of iqtest.d echo iqtest.f 1>&2 cat >iqtest.f <<'End of iqtest.f' PROGRAM IQTEST C TEST PROGRAM FOR IQPACK: DOUBLE PRECISION VERSION C DOUBLE PRECISION A,ALFA,ALPHA,B,BETA,FRQ,QFSUM,QFSX,T,WF,WTS INTEGER I,IER,IK,IND,IWF,KEY,KIND,LI,LO,LU,MLT,MTK,NDX,NIWF, 1 NT,NTMAX,NWF,NWTS DOUBLE PRECISION ZERO,HALF,ONE,TWO,THREE,FOUR,PI DIMENSION T(100),MLT(100),WTS(100),NDX(100),WF(10000),IWF(500) C PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TWO=2.0D0) PARAMETER (THREE=3.0D0,FOUR=4.0D0) PARAMETER (PI=3.14159265358979323846264338327950 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 C FUNCTIONS AND SUBROUTINES REFERENCED - CEGQF,CEGQFS,CEIQF, C CEIQFS,CGQF,CIQF,CIQFS,CLIQF,CLIQFS,EIQFS,F,COS, C SIN EXTERNAL F C LOGICAL UNIT NUMBERS: LO,LU=OUTPUT LI=INPUT (NOT USED) DATA LO/6/,LI/5/,LU/6/ C SIZE OF KNOTS ARRAY AND WORKSPACE IN DIMENSION STATEMENTS NTMAX=100 NWF=10000 NIWF=500 C GENERATE FEJER TYPE RULES - I.E. GAUSS CHEBYSHEV KNOTS AND C LEGENDRE WEIGHT FUNCTION C SET WEIGHT FUNCTION TYPE AND PARAMETERS KIND=1 A=ONE B=FOUR C ALPHA, BETA NOT USED IN LEGENDRE WEIGHT FUNCTION BUT SET ANYWAY ALPHA=HALF BETA=-HALF/TWO C C SET MULTIPLICITY OF KNOTS MTK=2 KEY=1 C NUMBER OF KNOTS. NT=5 C IF NT IS CHANGED IT MUST BE .LE. NTMAX C SET SIZE OF WEIGHTS ARRAY NWTS=MTK*NT FRQ=PI/TWO/NT C C TEST THE DRIVERS DO 110 IND=1,7 C C SET UP THE KNOTS FOR NT-POINT FEJER TYPE RULE C FOR IND=2,4,6 SCALE THE KNOTS TO GIVEN A,B C OTHERWISE - DEFAULT A,B - DON'T SCALE DO 10 I=1,NT T(I)=COS((TWO*I-1)*FRQ) IF(MOD(IND,2).EQ.0) THEN IF(KIND.EQ.5.OR.KIND.EQ.6.OR.KIND.EQ.8) THEN T(I)=T(I)+A ELSE T(I)=((B-A)*T(I)+(A+B))/TWO ENDIF ENDIF 10 CONTINUE C SET KNOT MULTIPLICITIES DO 20 I=1,NT MLT(I)=MTK 20 CONTINUE WRITE(LO,*)'----------------------------------------' GOTO(30,40,50,60,70,80,90),IND C 30 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CIQFS' CALL CIQFS(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA, 1 LU,NWF,WF,IER) WRITE(LO,*)' RETURN FROM CIQFS: ERROR INDICATOR IER= ',IER GOTO 110 C 40 CONTINUE C CIQF CALLS CIQFS, CGQF CALLS CGQFS C SO THIS TESTS BOTH WRITE(LO,*)' BEGINNING TEST OF CIQF,CIQFS,CGQF AND CGQFS' WRITE(LO,*)' WITH ALL CLASSICAL WEIGHT FUNCTIONS' A=-HALF B=TWO ALFA=-HALF DO 45 IK=1,8 KIND=IK BETA=TWO IF(KIND.EQ.8)BETA=-FOUR*FOUR WRITE(LO,*)' KNOTS AND WEIGHTS OF GAUSS QF' WRITE(LO,*)' COMPUTED BY CGQF(S)' CALL CGQF(NT,T,WTS,KIND,ALPHA,BETA,A,B,LU, 1 NWF,WF,NIWF,IWF,IER) WRITE(LO,*) WRITE(LO,*)' RETURN FROM CGQF(S): ERROR INDICATOR IER= ',IER C NOW COMPUTE THE WEIGHTS FOR THE SAME KNOTS BY CIQF WRITE(LO,*)' WEIGHTS OF GAUSS QF COMPUTED FROM THE' WRITE(LO,*)' KNOTS BY CIQF(S)' CALL CIQF(NT,T,MLT,NWTS,WTS,NDX,KEY,KIND,ALPHA,BETA,A,B, 1 LU,NWF,WF,IER) WRITE(LO,*) WRITE(LO,*)' RETURN FROM CIQF(S): ERROR INDICATOR IER= ',IER 45 CONTINUE GOTO 110 C 50 CONTINUE KIND=1 WRITE(LO,*)' BEGINNING TEST OF CEIQFS' QFSX=COS(-ONE)-COS(ONE) C QFSX IS EXACT VALUE OF INTEGRAL CALL CEIQFS(NT,T,MLT,KIND,ALPHA,BETA,F,QFSUM, 1 NWF,WF,NIWF,IWF,IER) WRITE(LO,*) WRITE(LO,*)' RETURN FROM CEIQFS: ERROR INDICATOR IER= ',IER WRITE(LO,*)' INTERGRAL OF SIN(X) ON [-1,1] BY FEJER TYPE RULE' WRITE(LO,*)' WITH ',NT,' POINTS OF MULTIPLICITY ',MTK WRITE(LO,*)' QUADRATURE FORMULA:',QFSUM WRITE(LO,*)' EXACT VALUE :',QFSX GOTO 110 C 60 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CEIQF' KIND=1 QFSX=COS(A)-COS(B) C QFSX IS EXACT VALUE OF INTEGRAL CALL CEIQF(NT,T,MLT,KIND,ALPHA,BETA,A,B,F,QFSUM, 1 NWF,WF,NIWF,IWF,IER) WRITE(LO,*) WRITE(LO,*)' RETURN FROM CEIQF: ERROR INDICATOR IER= ',IER WRITE(LO,*)' INTERGRAL OF SIN(X) FROM ',A,' TO ',B WRITE(LO,*)' BY FEJER TYPE RULE WITH ',NT,' POINTS' WRITE(LO,*)' OF MULTIPLICITY ',MTK WRITE(LO,*)' QUADRATURE FORMULA:',QFSUM WRITE(LO,*)' EXACT VALUE :',QFSX GOTO 110 C 70 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CLIQFS' CALL CLIQFS(NT,T,WTS,KIND,ALPHA,BETA,LU,NWF,WF,NIWF,IWF,IER) WRITE(LO,*) WRITE(LO,*)'RETURN FROM CLIQFS. IER= ',IER GOTO 110 C 80 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CLIQF AND EIQFS' CALL CLIQF(NT,T,WTS,KIND,ALPHA,BETA,A,B,LU,NWF,WF,NIWF,IWF,IER) WRITE(LO,*) WRITE(LO,*)' RETURN FROM CLIQF: ERROR INDICATOR IER= ',IER CALL EIQFS(NT,T,WTS,F,QFSUM,IER) WRITE(LO,*)' RETURN FROM EIQFS: ERROR INDICATOR IER= ',IER QFSX=COS(A)-COS(B) WRITE(LO,*)' INTERGRAL OF SIN(X) FROM ',A,' TO ',B WRITE(LO,*)' BY FEJER TYPE RULE WITH ',NT,' POINTS' WRITE(LO,*)' OF MULTIPLICITY ONE' WRITE(LO,*)' QUADRATURE FORMULA:',QFSUM WRITE(LO,*)' EXACT VALUE :',QFSX GOTO 110 90 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CEGQF' NT=12 C C EXPONENTIAL WEIGHT FUNCTION ON (A,B) CHOSEN C KIND=7 ALPHA=ONE BETA=ZERO A=ONE B=FOUR C C COMPUTE AND EVALUATE THE GAUSS QF CALL CEGQF(NT,KIND,ALPHA,BETA,A,B,F,QFSUM, 1 NWF,WF,NIWF,IWF,IER) WRITE(LO,*)' RETURN FROM CEGQF: ERROR INDICATOR IER= ',IER IF(IER.NE.0) GOTO 110 C QFSX=(B-A)*HALF*(COS(A)-COS(B))+SIN(B)+SIN(A) 1 -TWO*SIN((A+B)/TWO) WRITE(LO,*)' INTERGRAL OF X*SIN(X) FROM ',A,' TO ',B WRITE(LO,*)' BY GAUSS-EXPONENTIAL RULE WITH ',NT,' POINTS' WRITE(LO,*)' QUADRATURE FORMULA:',QFSUM WRITE(LO,*)' EXACT VALUE :',QFSX C 100 CONTINUE WRITE(LO,*)' BEGINNING TEST OF CEGQFS' NT=12 C C EXPONENTIAL WEIGHT FUNCTION, DEFAULT VALUES OF A,B C KIND=7 ALPHA=ONE BETA=ZERO C C COMPUTE AND EVALUATE THE GAUSS QF CALL CEGQFS(NT,KIND,ALPHA,BETA,F,QFSUM, 1 NWF,WF,NIWF,IWF,IER) WRITE(LO,*)' RETURN FROM CEGQFS: ERROR INDICATOR IER= ',IER IF(IER.NE.0) GOTO 110 C A=-ONE B=ONE QFSX=(B-A)*HALF*(COS(A)-COS(B))+SIN(B)+SIN(A) 1 -TWO*SIN((A+B)/TWO) WRITE(LO,*)' INTERGRAL OF X*SIN(X) FROM ',A,' TO ',B WRITE(LO,*)' BY GAUSS-EXPONENTIAL RULE WITH ',NT,' POINTS' WRITE(LO,*)' QUADRATURE FORMULA:',QFSUM WRITE(LO,*)' EXACT VALUE :',QFSX C 110 CONTINUE STOP END DOUBLE PRECISION FUNCTION F(X,I) C THIS FUNCTION GENERATES VALUES OF THE INTEGRAND C AND POSSIBLY ITS DERIVATIVES. SINCE IT APPEARS AS AN ACTUAL C PARAMETER IN THE ROUTINE WHICH USES IT, IT SHOULD BE DECLARED C IN AN EXTERNAL STATEMENT IN ANY CALLING PROGRAM DOUBLE PRECISION X INTEGER I,K,L C FUNCTIONS AND SUBROUTINES REFERENCED - COS,SIN,MOD C C RETURNS THE I-TH DERIVATIVE OF SIN(T) AT T=X C (ZERO-TH DERIVATIVE IS THE FUNCTION VALUE) K=MOD(I,4) L=MOD(K,2) IF(L.EQ.0)F=SIN(X) IF(L.EQ.1)F=COS(X) IF(K.GT.1)F=-F RETURN END End of iqtest.f