C ALGORITHM 672, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 15, NO. 2, PP. 137-143. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration driver C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C The following templates demonstrate the use of the procedure EXTEND C to generate sequences of extended quadrature rules for various weight C functions given the definitions of the orthogonal polynomial 3-term C recurrence relations. ICASE selects the C required sequence as follows: C ICASE = 1 3-point Gauss-Legendre in [-1,1] C ICASE = 2 2-point Gauss-Lobatto in [-1,1] C ICASE = 3 6-point Radau in [-1,1] C ICASE = 4 2-point Gauss-Laguerre in [0,infinity) C ICASE = 5 3-point Gauss-Hermite in (-infinity,infinity) C ICASE = 6 3-point Gauss-Jacobi in [0,1] C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C FORTRAN-77 code C Unless indicated otherwise the type of each variable is implied w%C by the def C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (LDA=257,LDB=2*LDA+1,IDIGIT=8,NEXP=38,NTEST=4) PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0) DOUBLE PRECISION T(0:LDA),EXT(0:LDA) DOUBLE PRECISION QR(LDA),QI(LDA) DOUBLE PRECISION PNODES(LDA),WT(LDA) DOUBLE PRECISION ERR(LDA),TEST(0:NTEST) DOUBLE PRECISION WORKA(LDA,LDA),WORKB(LDB,3) INTEGER IWORK(LDA) LOGICAL SYMMET,START EXTERNAL RECURA,RECURB,RECURC,RECURD C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Get data from terminal C C Select demonstration 5 WRITE(*,10) 10 FORMAT(' Case ?:') READ(*,*)ICASE ICASE=MAX(1,ICASE) C Select number of iterative extensions to be performed WRITE(*,20) 20 FORMAT(' No. of rules?:') READ(*,*)NSEQ NSEQ=MAX(1,NSEQ) I=1 50 GOTO(100,200,300,400,500,600),ICASE C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 1 - extension of 3-point Gauss-Legendre rule C Recurrence defined by RECURA 100 IF(I.EQ.1) THEN C Generate 3-point Gauss initially from zero point rule, C i.e. no pre-assigned nodes, symmetry exploited. WRITE(*,*) 'Gauss-Legendre 3-point extension' N=0 M=3 M0=0 T(0)=ONE SYMMET=.TRUE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 120 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 120 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 2 - extension of 2-point Lobatto rule C Recurrence defined by RECURA 200 IF(I.EQ.1) THEN C Add one node, Pre-assign -1.0 and 1.0, symmetry exploited. WRITE(*,*) 'Lobatto 2-point extension' N=2 M=1 PNODES(1)=ONE PNODES(2)=-ONE SYMMET=.TRUE. START=.TRUE. ELSE C Add N-1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N-1 START=.FALSE. END IF H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 220 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 220 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 3 - extension of 6-point Radau rule C Recurrence defined by RECURA 300 IF(I.EQ.1) THEN C Add five nodes. Pre-assign -1.0, no symmetry. WRITE(*,*) 'Radau 6-point extension' N=1 M=5 PNODES(1)=-ONE SYMMET=.FALSE. START=.TRUE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 320 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 320 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 4 - extension of 2-point Gauss-Laguerre C Recurrence defined by RECURB 400 IF(I.EQ.1) THEN C Generate 2-point rule initially from zero point rule, C i.e. no pre-assigned nodes, no symmetry. WRITE(*,*) 'Gauss-Laguerre 2-point extension' N=0 M=2 M0=0 T(0)=ONE SYMMET=.FALSE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=ONE IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURB,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 420 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURB,TEST(K),IERR) 420 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 5 - extension of 3-point Gauss-Hermite rule C Recurrence defined by RECURC 500 IF(I.EQ.1) THEN C Generate 3-point rule initially from zero point rule, C i.e. no pre-assigned nodes, symmetry exploited. WRITE(*,*) 'Gauss-Hermite 3-point extension' M=3 N=0 M0=0 T(0)=ONE SYMMET=.TRUE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension C Zero moment integral = sqrt(pi) H0=TWO*SQRT(ATAN(ONE)) IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURC,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 520 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURC,TEST(K),IERR) 520 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Demonstration 6 - extension of 3-point Gauss-Jacobi C for weight sqrt(x) in [0,1] C Recurrence defined by RECURD 600 IF(I.EQ.1) THEN C Generate 3-point rule initially from zero point rule, C i.e. no pre-assigned nodes, no symmetry. WRITE(*,*) 'Gauss-Jacobi 3-point extension' M=3 N=0 M0=0 T(0)=ONE SYMMET=.FALSE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=TWO/THREE IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURD,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 620 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURD,TEST(K),IERR) 620 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C Display results 2000 WRITE(*,*)'Iteration',I WRITE(*,*)'Coefficients of expansion whose roots ', * 'are the new nodes:' WRITE(*,3000)(EXT(J),J,J=0,M) C WRITE(*,*)'New nodes' WRITE(*,3100)(QR(K),QI(K),IWORK(K),ERR(K),K=1,NODES) C WRITE(*,*)'New full extended expansion' WRITE(*,3010)(T(J-M0),J,J=M0,N) C WRITE(*,3200)I,N,IFLAG,NODES IF(IFLAG.NE.0.AND.IFLAG.NE.6) THEN WRITE(*,*)'Terminated prematurely - see IFLAG' GOTO 5 END IF C Print rule (positive nodes only if symmetry present) IF(SYMMET) THEN NUM=(N+1)/2 ELSE NUM=N END IF WRITE(*,3300)(J,PNODES(J),WT(J),J=1,NUM) C Display test results DO 2010 K=0,MIN(NTEST,IDEG/2) WRITE(*,3400)K,TEST(K) 2010 CONTINUE IF(IFLAG.EQ.6) THEN WRITE(*,*)'Rule test is unsatisfactory' GOTO 5 END IF I=I+1 IF(I.LE.NSEQ) GOTO 50 GOTO 5 C 3000 FORMAT(D25.16,'*P(',I3,',X)') 3010 FORMAT(D25.16,'*P(',I3,',X)/HI') 3100 FORMAT(21X,'REAL',16X,'IMAGINARY',1X,'FLAG',7X,'ERR'/, * (2D25.16,I5,D10.1)) 3200 FORMAT(' Complete extended rule: STEP=',I2,' POINTS=',I3, * ' IFLAG=',I1,' NODES ADDED=',I3) 3300 FORMAT(2X,'No.',21X,'NODE',19X,'WEIGHT',/,(I5,2D25.16)) 3400 FORMAT(' TEST(',I2,')=',D25.16) END C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C User defined subroutines - Recurrence relations C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- SUBROUTINE RECURA(K,C,D,E) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C This is an example of a user supplied subroutine to define the C orthogonal polynomials. C C CALL RECUR(K,C,D,E) gives the coefficients C,D and E such that, C C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameter: C K = Index for relation C Output parameters: C C,D,E = Parameters in the recurrence relation C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0) C Legendre recurrence for [-1,1] C Covers Gauss, Lobatto and Radau F=FLOAT(K+1) C=FLOAT(2*K+1)/F D=ZERO E=-FLOAT(K)/F RETURN END SUBROUTINE RECURB(K,C,D,E) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0D0) C Laguerre recurrence F=FLOAT(K+1) C=-ONE/F D=FLOAT(2*K+1)/F E=-FLOAT(K)/F RETURN END SUBROUTINE RECURC(K,C,D,E) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,ONE=1.0D0,TWO=2.0D0) C Hermite recurrence C=ONE D=ZERO E=-FLOAT(K)/TWO RETURN END SUBROUTINE RECURD(K,C,D,E) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0D0,TWO=2.0D0,ONEP5=1.5D0) C Jacobi polynomials in [0,1] C Weight function is (1-x)**(p-q)*x**(q-1) C C Case for weight sqrt(x) i.e. p=3/2 and q=p P=ONEP5 Q=P FK=FLOAT(K) F2K=FK*TWO B=F2K+P BP1=B+ONE X3=(B-TWO)*(B-ONE)*B A1=(B-ONE)*BP1*X3 A2=-X3*(F2K*(FK+P)+Q*(P-ONE)) A3=X3*BP1*(B-ONE) A4=FK*(FK+Q-ONE)*(FK+P-ONE)*(FK+P-Q)*BP1 C=A3/A1 D=A2/A1 E=-A4/A1 RETURN END C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Main algorithm begins here with structure: C -> ASSIGN C C -> GENER ----> EPROD C C EXTEND --------> SOLVE ----> NEWTON/BAIR C v C v -> RSORT C v C CHECK -> WEIGHT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- SUBROUTINE EXTEND(N,M,M0,T,RECUR,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QRNODE,QINODE,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION T(0:*),EXT(0:*) DOUBLE PRECISION PNODES(*),WT(*),H0 DOUBLE PRECISION QRNODE(*),QINODE(*),ERR(*) DOUBLE PRECISION WORKA(0:LDA-1,0:*),WORKB(0:LDB-1,*) INTEGER M0,N,M,IWORK(*),LDA,LDB,NODES,IFLAG,NEXP,IDIGIT LOGICAL SYMMET,START EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Please address queries or comments to: C C T.N.L. Patterson C Department of Applied Mathematics & Theoretical Physics C The Queen's University of Belfast C Belfast, BT9 1NN C N. Ireland C C Tel: International +44 232 245133 Ext. 3792. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the N+M node quadrature rule composed of N pre-assigned nod C together with M nodes chosen optimally to achieve algebraic degree of C precision of at least N+2*M-1. C C The orthogonal system of polynomials associated with the quadrature C weight is defined generally by the recurrence relation specified in th C user supplied subroutine RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C N = Number of pre-assigned nodes (and upper limit to the expansion C Note that if successful this is reset to N+M on completion C (the appropriate value for iterative use). C M = Number of nodes to be optimally added. C M0 = Lower limit to the expansion of T. This is ignored if START C is .TRUE. Note that if successful this is reset to M C on completion (the appropriate value for iterative use). C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M0 to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. Element T(I-M0) holds the C value of TI. C C Note that T is either, C (1) provided explicitly, C (2) generated automatically from the N pre-assigned nodes C given in PNODES(*) (if START is .TRUE.) C or, C (3) generated from a previous call to the subroutine. C This array should be declared to have at least C max(N-M0+1,M+1) elements in the calling program. C C The service subroutine TRANSF can be used to transform C the expansion to the required input form if desired C with the parameter IFLAG set to 1. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C SYMMET = .FALSE. if no advantage is to be taken of symmetry, if any, C about x=0 in the interval of integration and the C orthogonality weight function. Note that if symmetry in C fact does exist setting this parameter to zero will still C produce correct results - only efficiency is effected. C = .TRUE. if the extended rule computations should C exploit symmetry about x=0 in the interval of C integration and the orthogonality weight function. C This reduces the size of the system of linear equations C determining EXT by a factor of about 2 (see WORKA). If C symmetry does not in fact exist erroneous results will be C produced. C START = .TRUE. then the polynomial T is generated to have C the pre-assigned nodes (PNODES) as its roots. C = .FALSE. then the supplied values of the coefficients C of T are used directly (see above). C PNODES = Array holding the pre-assigned nodes. This array should C be declared to have at least N+M elements in the calling prog C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. C NEXP = Largest negative decimal exponent supported on the C computer. (Positive number - typical value 38 for VAX/VMS). C Weights less than approximately 10**(-NEXP) are set to zero C when the Christoffel-Darboux identity is used (N=M). C This may be set to INT(LOG10(X1MACH(2))) where X is set to C correspond to the appropriate precision in the PORT library. C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C IWORK = Integer working array which should be declared in the C calling program to have at least max(M,N) elements. C On return IWORK provides information on the convergence C of the nodes. See output parameters. C WORKA = Real working matrix which should be declared in the calling C program to have dimension at least max(M+1,N) C by max(M+1,N+1). If SYMMET=.TRUE. (see above) the C dimension can be reduced to max(M/2+1,N) C by max(M/2+1,N+1). C LDA = Number of elements in the leading dimension of WORKA C declared in the calling program. C WORKB = Real working matrix which should be declared in the calling C program to have dimension at least 2*M+1 by 3. C LDB = Number of elements in the leading dimension of WORKB C declared in the calling program C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C PNODES = Ordered array holding the N+M nodes of the extended C quadrature rule made up from the original pre-assigned C nodes and the new optimally extended nodes. These values can C be used in subsequent iterative use of the subroutine. C WT = Array holding the values of the quadrature weights for C the extended rule associated with the nodes held in PNODES. C This array should be declared to have at least N+M elements C in the calling program. C T = Array holding the coefficients TI of the new orthogonal C expansion whose roots are the nodes of the extended quadratur C (that is, the pre-assigned nodes plus the extended nodes) and C is expressed as: C SUM (I=M to N+M) (TI/HI)*P(I,X) C T(I-M) holds the value of TI. C (For definitions see description of input argument T). C This polynomial can be used as input for further extensions. C C The service subroutine TRANSF can be used to remove the C moment factors from the expansion if desired with the C parameter IFLAG set to 0. C M0 = Lower limit defining the new orthogonal expansion T. C (Set to M). C N = Upper limit defining the new orthogonal expansion T. C (Set to the original value of N+M). C NODES = Number of extended nodes found. Normally equals M but see IFL C QRNODE = Array holding the real parts of the extended nodes (1,..,NODE C This array should be declared to have at least M elements C in the calling program. C QINODE = Array holding the imaginary parts of the extended C nodes (1,..,NODES). (Hopefully these values are zero!). C This array should be declared to have at least M elements C in the calling program. C ERR = Array holding a measure of the relative error in the C nodes. This may be inspected if the convergence C error flag has been raised (IFLAG=3) to decide if the nodes C in question are acceptable. (ERR(*) actually gives the mean C last correction to the quadratic factor in the generalised C Bairstow root finder (see BAIR). This should declared in C the calling program to have at least M elements. C EXT = Array holding the coefficients of the polynomial whose C roots are the extended nodes (QRNODES(*),QINODES(*)) and C expressed as: C EXT = SUM (I=0 to M) EXT(I)*P(I,X) C This array should be declared to have at least M+1 elements C in the calling program. C IWORK = Node convergence flags. Elements 1 to NODES give information C on the convergence of the roots of the polynomial EXT C corresponding to each extended node. C Element I = 0 Convergence of I th root satisfactory C Element I = 1 Convergence of I th root unsatisfactory C IFLAG = 0, No error detected C = 1, The linear system of equations defining the polynomial C whose roots are the extended nodes became singular or C very ill-conditioned. (FATAL). C = 2, The linear system of equations used to generate the C polynomial T when START is .TRUE. became singular C or very ill-conditioned. (FATAL). C = 3, Poor convergence has been detected in the calculation C of the roots of EXT (see above) corresponding to the new C nodes or all nodes have not been found (M not equal C to NODES). See also ERR(*) below. C = 4, Possible imaginary nodes detected. C = 5, Value of N and M incompatible for SYMMET=.TRUE. C Both cannot be odd. (FATAL) C = 6, Test of new quadrature rule has failed. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are included and renamed GEFA77 and GESL7 C These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are renamed DSCAL7, IDAMX7, DAXPY7 and DDOT7. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Changing the precision. C C This is accomplished as follows: C (1) Amend the TYPE statements, C (2) Select an appropriate value for the NEXP argument to EXTEND. C C NOTE: C (a) All constants used are specified in PARAMETER statements at the st C of each subprogram and, C (b) Generic names are used for all function calls. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- INTEGER S C IFLAG=0 NODES=0 IDEG=N+2*M-1 C Look for incompatible values of N and M IF(SYMMET) THEN C Both N and M cannot be odd IF(MOD(N,2).EQ.1.AND.MOD(M,2).EQ.1) THEN IFLAG=5 RETURN END IF END IF C Generate if required the initial T polynomial corresponding to C prescribed pre-assigned nodes IF(START .AND. N.NE.0) THEN CALL ASSIGN(N,PNODES,IWORK,WORKA,LDA,RECUR,T,IERR) M0=0 IF(IERR.NE.0) THEN IFLAG=2 RETURN END IF END IF NLAST=N C Generate extended expansion coefficients and overwrite T CALL GENER(T,M0,N,M,RECUR,SYMMET,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IERR) IF(IERR.NE.0) THEN IFLAG=1 RETURN END IF C Find extended nodes as roots of EXT(*) CALL SOLVE(EXT,M,SYMMET,RECUR,IDIGIT,QRNODE,QINODE, * NODES,ERR,IWORK,WORKB,LDB,IERR) IF(IERR.NE.0) IFLAG=IERR+2 IF(IFLAG.NE.0) RETURN C Accumulate nodes for extended rule DO 10 I=1,M PNODES(NLAST+I)=QRNODE(I) 10 CONTINUE C Re-order CALL RSORT(PNODES,N,1) C Compute weights (only for positive nodes if symmetric) IF(SYMMET) THEN NUM=(N+1)/2 ELSE NUM=N END IF DO 20 I=1,NUM CALL WEIGHT(T,M0,N,PNODES(I),RECUR,H0,NEXP,WT(I)) IF(SYMMET) THEN WT(N-I+1)=WT(I) END IF 20 CONTINUE C Test the new rule DO 30 K=0,MIN(4,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECUR,TEST,IERR) IF(IERR.EQ.1) THEN IFLAG=6 RETURN END IF 30 CONTINUE RETURN END SUBROUTINE CHECK(N,QNODE,WT,K,H0,RECUR,TEST,IERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION QNODE(*),WT(*),H0,TEST INTEGER N,K,IERR EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Carry out test of the given quadrature rule by computing the C appropriate integral of, C W(X)*P(K,X)*P(K,X) C over the region associated with the weight function W(X) and the C orthogonal polynomials P(K,X). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Number of nodes in the quadrature rule. C QNODE = Array holding the N nodes. C WT = Array holding the N weights. C K = Index of the orthogonal polynomial whose weighted square C is to be integrated. C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. Note that C P(0,X) is arbitrarily taken to be 1.0. C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C C Output parameters: C TEST = Approximate value of the test integral normalised to C unity. Thus, ABS(TEST-1) gives a measure of the C quality of the calculated rule. C IERR = 0, OK. C = 1, Rule quality unsatisfactory C = 2, Invalid values for input arguments C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,ONE=1.0D0,TOL=0.0000001D0) IERR=0 IF(K.LT.0 .OR. N.LT.1 .OR. H0.LE.ZERO) THEN IERR=2 RETURN END IF TEST=ZERO DO 30 I=1,N P1=ONE IF(K.EQ.0) GOTO 20 P0=ZERO X=QNODE(I) C Calculate integrand DO 10 J=0,K-1 CALL RECUR(J,CJ,DJ,EJ) P=(CJ*X+DJ)*P1+EJ*P0 P0=P1 P1=P 10 CONTINUE 20 TEST=TEST+P1*P1*WT(I) 30 CONTINUE TEST=TEST/H0 IF(K.EQ.0) RETURN C Calculate exact value CALL RECUR(0,P,P0,P1) DO 70 J=1,K CALL RECUR(J,CJ,DJ,EJ) P=-P*EJ 70 CONTINUE C Normalise result to unity TEST=TEST*CJ/P C Test for rule quality IF(ABS(TEST-ONE).GT.TOL) IERR=1 RETURN END SUBROUTINE ASSIGN(N,PNODES,IWORK,WORK,LDW,RECUR,T,IERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION PNODES(*),WORK(0:LDW-1,0:*),T(0:*) INTEGER N,LDW,IERR,IWORK(*) EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Generate the initial polynomial T whose roots are the required C pre-assigned nodes C C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C N = Number of pre-assigned nodes to be used to generate T. C PNODES = Array holding N pre-assigned nodes to be be used to C generate T. C IWORK = Integer working array which should be declared in the C calling program to have at least N elements. C WORK = Real working matrix which should be declared in C the calling program to have dimension at least N by C N+1. C LDW = Number of elements in the leading dimension of WORK C declared in the calling program C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C T = Array holding the coefficients of the polynomial whose C roots define the pre-assigned nodes of the quadrature C rule and expressed as: C H0* SUM (I=0 to N) T(I)/HI*P(I,X) C T(I) holds the value of TI. C This array should be declared to have at least N+1 elements C in the calling program. C IERR = 0, No error detected C = 1, The linear system of equations used to generate the C polynomial T became singular or very ill-conditioned. (FAT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C External library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are used and renamed GEFA77 and GESL77. C Note -- These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are not renamed. C (Quadruple precision versions used for this subprogram) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,ONE=1.0D0) IERR=0 C Set up the linear system of equations DO 20 J=1,N X=PNODES(J) P0=ZERO P1=ONE P=P1 DO 10 K=0,N WORK(J-1,K)=P CALL RECUR(K,C0,D0,E0) P=(C0*X+D0)*P1+E0*P0 P0=P1 P1=P 10 CONTINUE 20 CONTINUE C Solve linear system CALL GEFA77(WORK,LDW,N,IWORK,INFO) IF(INFO.NE.0) THEN IERR=1 RETURN END IF CALL GESL77(WORK,LDW,N,IWORK,WORK(0,N),0) DO 30 J=0,N-1 T(J)=-WORK(J,N) 30 CONTINUE T(N)=ONE C Weight with moments CALL TRANSF(T,0,N,RECUR,1) RETURN END SUBROUTINE GENER(T,M0,N,M,RECUR,SYMMET,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION WORKA(0:LDA-1,0:*),WORKB(0:LDB-1,*) DOUBLE PRECISION T(0:*),EXT(0:*) INTEGER M0,N,M,IWORK(*),LDA,LDB,IFLAG LOGICAL SYMMET C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Given N pre-assigned quadrature nodes defined as the roots of the C polynomial expansion, C SUM (I=M0 to N) (TI/HI)*P(I,X) C calculate the polynomial expansion, C SUM (I=0 to M) SI*P(I,X) C whose roots are the M optimal nodes and new expansion C SUM (I=M to N+M) (RI/HI)*P(I,X) C whose roots are to the N+M nodes of the full extended quadrature rule. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M0 to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. T(I-M0) holds the C value of TI. This array should be declared to have at least C max(N-M0+1,M+1) elements in the calling program. C M0 = Lower limit to the expansion of T. C N = Upper limit to expansion of T. C M = Number of nodes to be optimally added. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C SYMMET= .FALSE. if no advantage is to be taken of symmetry, if any, C about x=0 in the interval of integration and the C orthogonality weight function. Note that if symmetry in C fact does exist setting this parameter to zero will still C produce correct results - only efficiency is effected. C = .TRUE. if the extended rule computations should C exploit symmetry about x=0 in the interval of C integration and the orthogonality weight function. C This reduces the size of the system of linear equations C determining EXT by a factor of about 2 (see WORKA). If C symmetry does not in fact exist erroneous results will be C produced. C IWORK = Integer working array which should be declared in the C calling program to have at least M elements. C WORKA = Real working matrix which should be declared in the calling C program to have dimension at least M+1 by max(M+1,N+1). C If SYMMET=.TRUE. (see above) the dimension can be reduced to C M/2+1 by max(M/2+1,N/2+1). C LDA = Number of elements in the leading dimension of WORKA C declared in the calling program C WORKB = Real working matrix which should be declared in the calling C program to have at dimension at least 2*M+1 by 3. C LDB = Number of elements in the leading dimension of WORKB C declared in the calling program. C C Output parameters: C T = Array holding the coefficients of the new orthogonal C expansion whose roots are the nodes of the extended quadratur C (that is the pre-assigned nodes plus the extended nodes). C It is expressed as: C SUM (I=M to N+M) (TI/HI)*P(I,X) C where N and M have their original values. T(I-M) holds C the value of TI. See input argument of T for definitions. C M0,N = Lower and upper limits defining the new orthogonal expansion C EXT = Array holding the coefficients of the polynomial whose C roots are the new extended nodes and expressed as: C EXT = SUM (I=0 to M) EXT(I)*P(I,X) C IFLAG = 0, No error detected C = 1, The linear system of equations defining the polynomial C whose roots are the extended nodes became singular or C very ill-conditioned. (FATAL). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C External library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are used and renamed GEFA77 and GESL77. C Note -- These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are not renamed. C (Quadruple precision versions used for this subprogram) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- LOGICAL NEVEN,MSODD,MISS EXTERNAL RECUR INTEGER S PARAMETER (ZERO=0.0D0,ONE=1.0D0) IFLAG=0 C Look for trivial case IF(N.EQ.0) THEN DO 10 I=0,M-1 EXT(I)=ZERO 10 CONTINUE EXT(M)=ONE T(0)=ONE N=M M0=M RETURN END IF C General case NEVEN=MOD(N,2).EQ.0 NM=N+M C Form matrix 20 DO 60 S=0,M MSODD=MOD(M+S,2).EQ.1 IF(NEVEN.AND.MSODD.AND.SYMMET) GOTO 60 DO 50 J=0,S CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IF(MOD(N+S+J,2).EQ.1.AND.SYMMET) GOTO 50 IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 40 DO 30 I=IBOT,ITOP SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 30 CONTINUE 40 IF(.NOT.SYMMET) THEN WORKA(S,J)=SUM WORKA(J,S)=SUM GOTO 50 END IF IF(NEVEN) THEN WORKA(S/2,J/2)=SUM WORKA(J/2,S/2)=SUM ELSE IF(MSODD) THEN WORKA(S/2,J/2)=SUM ELSE WORKA(J/2,S/2)=SUM END IF END IF 50 CONTINUE 60 CONTINUE NEQ=M IF(SYMMET) NEQ=M/2 C Solve for expansion coefficients CALL GEFA77(WORKA,LDA,NEQ,IWORK,INFO) IF(INFO.NE.0) THEN IFLAG=1 RETURN END IF CALL GESL77(WORKA,LDA,NEQ,IWORK,WORKA(0,NEQ),0) C Store expansion coefficients DO 70 J=0,NEQ-1 EXT(J)=-WORKA(J,NEQ) 70 CONTINUE EXT(NEQ)=ONE C Calculate new T polynomial IF(SYMMET) GOTO 160 C C Non-symmetric case DO 140 S=M,NM IF(S.EQ.M) GOTO 120 DO 110 J=0,M CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 100 DO 90 I=IBOT,ITOP IR=I-IREF SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 90 CONTINUE 100 WORKA(M,J)=SUM 110 CONTINUE 120 SUM=ZERO DO 130 I=0,M SUM=SUM+EXT(I)*WORKA(M,I) 130 CONTINUE WORKA(M-1,S-M)=SUM 140 CONTINUE C Overwrite old values of T DO 150 I=0,N T(I)=WORKA(M-1,I) 150 CONTINUE GOTO 250 C C Symmetric case 160 DO 210 S=M,NM IF(MOD(N+M+S,2).EQ.1) GOTO 210 DO 190 J=0,M CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IF(MOD(N+S+J,2).EQ.1) GOTO 190 IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 180 DO 170 I=IBOT,ITOP IR=I-IREF SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 170 CONTINUE 180 WORKA(NEQ,J/2)=SUM 190 CONTINUE SUM=ZERO DO 200 I=0,NEQ SUM=SUM+EXT(I)*WORKA(NEQ,I) 200 CONTINUE WORKA(NEQ-1,(S-M)/2)=SUM 210 CONTINUE C Overwrite old values of T in full unsymmetric form IC=N/2 MISS=.TRUE. DO 220 J=N,0,-1 MISS=.NOT.MISS IF(MISS) THEN T(J)=ZERO ELSE T(J)=WORKA(NEQ-1,IC) IC=IC-1 END IF 220 CONTINUE C Convert EXT to full unsymmetric form WORKB(M,1)=ONE IC=NEQ-1 MISS=.FALSE. DO 230 J=M-1,0,-1 MISS=.NOT.MISS IF(MISS) THEN WORKB(J,1)=ZERO ELSE WORKB(J,1)=EXT(IC) IC=IC-1 END IF 230 CONTINUE DO 240 J=0,M EXT(J)=WORKB(J,1) 240 CONTINUE C Scale new T polynomial 250 PMAX=ZERO DO 260 I=0,N PMAX=MAX(PMAX,ABS(T(I))) 260 CONTINUE DO 270 I=0,N T(I)=T(I)/PMAX 270 CONTINUE N=NM M0=M RETURN END SUBROUTINE SOLVE(EXT,M,SYMMET,RECUR,IDIGIT,QRNODE,QINODE, * NODES,ERR,ICHECK,WORK,LDW,IERR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION EXT(0:*),WORK(0:LDW-1,*),ERR(*) DOUBLE PRECISION QRNODE(*),QINODE(*) INTEGER M,NODES,LDW,IERR,ICHECK(*),IDIGIT LOGICAL SYMMET,RESET EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculate the roots of the orthogonal polynomial expansion C expressed as, C SUM (I=0 to M) EXT(I)*P(I,X) C where the array EXT holds the appropriate coefficients. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C EXT = Array holding the coefficients of the polynomial whose C roots are required (nodes of the quadrature rule) C and expressed as: C SUM (I=0 to M) EXT(I)*P(I,X) C The recurrence relation for the orthogonal polynomials C P(I,X) is defined by the subroutine RECUR. C This array should be declared to have at least M+1 elements C in the calling program. C M = Upper limit to expansion EXT (polynomial degree). C SYMMET = .FALSE. if no advantage can be taken of symmetry C about x=0 in the interval of integration and the C orthogonality weight function. C = .TRUE. if symmetry exists about x=0 in the interval of C integration and the orthogonality weight function. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C WORK = Real working matrix which should be declared in the calling C program to have dimension at least M+1 by 2. C LDW = Number of elements in the leading dimension of WORK C declared in the calling program C C Output parameters: C C QRNODE = Array holding the real parts of the roots fo EXT (1,..,NODES) C QINODE = Array holding the imaginary parts of the roots of EXT (1,..,N C (hopefully these values are zero!). C NODES = Number of extended nodes found. Normally equals M but see IER C ICHECK = Root convergence flags. Elements 1 to NODES give information C on the convergence of the roots of the polynomial EXT. C Element I = 0 Convergence of I th root satisfactory C Element I = 1 Convergence of I th root unsatisfactory C This array should be declared to have at least M elements C in the calling program. C IERR = 0, No error detected C = 1, Possible imaginary nodes detected. C = 2, Poor convergence has been detected in the calculation C of the roots of EXT (see above) or all roots have not C been found (M not equal to NODES). See also ERR(*) below. C ERR = Array holding a measure of the relative error in the C roots. This may be inspected if the convergence C error flag has been raised (IERR=2) to decide if the roots C in question are acceptable. (ERR(*) actually gives the mean C last correction to the quadratic factor in the generalised C Bairstow root finder (see BAIR). This should declared in C the calling program to have at least M elements. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0, * TWO=2.0D0,VRT1=0.0000001D0) NODES=0 IERR=0 C If M is odd find and remove initial real root using NEWTON iteration C Set WORK(*,1) to polynomial to be processed IF(MOD(M,2).EQ.1) THEN ZR1=VRT1 CALL NEWTON(EXT,M,ZR1,RECUR,IDIGIT,WORK(0,1),ERRVAL,IFAIL) NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR1 QINODE(NODES)=ZERO NROOT=M-1 ELSE DO 10 I=0,M WORK(I,1)=EXT(I) 10 CONTINUE NROOT=M END IF IF(NROOT.EQ.0) GOTO 50 C Find remaining root pairs C Calculate seed approximation for quadratic factor CALL RECUR(0,C0,D0,E0) CALL RECUR(1,C1,D1,E1) RT1=VRT1 RT2=ZERO IF(SYMMET) RT2=-RT1 P1A=C0*RT1+D0 P1B=C0*RT2+D0 P2A=(C1*RT1+D1)*P1A+E1 P2B=(C1*RT2+D1)*P1B+E1 DET=C0*(RT1-RT2) SA1=(P2A-P2B)/DET SA0=(P1A*P2B-P1B*P2A)/DET RESET=.TRUE. C Alternate approximation which introduces small complex component RT1=VRT1 RT2=VRT1 SFA1=(C0*D1+D0*C1)/C0+TWO*C1*RT1 SFA0=D0*D1+E1-D0*SFA1-C0*C1*(RT1*RT1+RT2*RT2) C IP1 points to current deflated polynomial IP1=1 IP2=2 20 IF(RESET) THEN A2=ONE A1=SA1 A0=SA0 RESET=.FALSE. END IF CALL BAIR(NROOT,WORK(0,IP1),WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) IF(IFAIL.NE.0) THEN C Try again with complex components introduced A2=ONE A1=SFA1 A0=SFA0 RESET=.TRUE. CALL BAIR(NROOT,WORK(0,IP1),WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) END IF C Apply Bairstow to full expansion to avoid error accumulation CALL BAIR(M,EXT,WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) C Tidy up the quotient polynomial CALL QFACT(NROOT,WORK(0,IP1),WORK(0,IP2),RECUR,A1,A0, * ZR1,ZR1,ZR1,ZR1,ZR1,ZR1) CALL ROOTS(A0,A1,A2,ZR1,ZI1,ZR2,ZI2,RECUR,INFO) C Record node information NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR1 QINODE(NODES)=ZI1 NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR2 QINODE(NODES)=ZI2 NROOT=NROOT-2 C Make the deflated polynomial current I=IP1 IP1=IP2 IP2=I IF(NROOT.GT.0) THEN C Scale the deflated polynomial PMAX=ZERO DO 30 I=0,NROOT PMAX=MAX(PMAX,ABS(WORK(I,IP1))) 30 CONTINUE DO 40 I=0,NROOT WORK(I,IP1)=WORK(I,IP1)/PMAX 40 CONTINUE GOTO 20 END IF C Calculation complete - Check for difficulties C Look for poor convergence 50 I=0 DO 60 J=1,NODES I=I+ICHECK(J) 60 CONTINUE IF(NODES.NE.M.OR.I.NE.0) THEN IERR=1 RETURN END IF C Look for possible imaginary nodes DO 70 J=1,NODES IF(QINODE(J).NE.ZERO) THEN IERR=2 RETURN END IF 70 CONTINUE RETURN END SUBROUTINE RSORT(A,N,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(*) INTEGER N,IFLAG C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C Carries out a simple ripple sort of A(*). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C A = Array holding the numbers to be sorted C N = Number of elements to be sorted C IFLAG = 0 for ascending sort C = 1 for descending sort C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- LOGICAL DONE,ASCEND ASCEND=IFLAG.EQ.0 C Begin scans DO 30 J=N-1,1,-1 DONE=.TRUE. DO 20 K=1,J IF(ASCEND) THEN K1=K K2=K+1 ELSE K1=K+1 K2=K END IF IF(A(K1).GT.A(K2)) THEN C Exchange elements VAL=A(K1) A(K1)=A(K2) A(K2)=VAL DONE=.FALSE. END IF 20 CONTINUE IF(DONE) RETURN 30 CONTINUE RETURN END SUBROUTINE WEIGHT(T,M,N,XNODE,RECUR,H0,NEXP,WT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION T(0:*),XNODE,H0,WT INTEGER M,N,NEXP EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the quadrature weight associated with the node XNODE in the C rule whose nodes are defined by the roots of polynomial T. C C The weight is calculated by dividing T by (X-XNODE) to give, C C S(X) = T(X)/(X-XNODE) = SUM (0 to N-1) G(I)*P(I,X). C C S(X) is then divided by (X-XNODE) to give remainder R. C C The weight is finally given by H0*G(0)/R. If N=M the C Christoffel-Darboux identity result is used to reduce extreme C cancellation effects at high degree. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. T(I-M) holds the C value of TI. This array should be declared to have at least C N-M+1 elements in the calling program. C M = Lower limit to the expansion of T C N = Upper limit to expansion of T C XNODE = Node whose weight is required C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. Note that C P(0,X) is arbitrarily taken to be 1.0 C NEXP = Largest negative decimal exponent supported on the C computer. (Positive number - typical value 38). C Weights less than approximately 10**(-NEXP) are set to zero C when the Christoffel-Darboux identity is used (N=M). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C WT = Weight associated with XNODE. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,ONE=1.0D0,TEN=10.0D0) C Check for special case IF(M.EQ.N) THEN C Use Christoffel-Darboux result BK1=ZERO BK2=ONE DK1=ZERO DK2=ZERO ISCALE=0 CALL RECUR(0,H,D0,E0) DO 20 K=0,N-1 CALL RECUR(K,CK,DK,EK) IF(K.NE.0) H=-EK*H BB=(CK*XNODE+DK)*BK2+EK*BK1 DD=(CK*XNODE+DK)*DK2+EK*DK1+CK*BK2 BK1=BK2 BK2=BB DK1=DK2 DK2=DD IF(BK2.NE.ZERO) THEN J=LOG10(ABS(BK2)) IF(ABS(J).GT.2) THEN C Scale to control overflow/underflow ISCALE=ISCALE-2*J SCALE=TEN**J BK2=BK2/SCALE BK1=BK1/SCALE DK1=DK1/SCALE DK2=DK2/SCALE END IF END IF IF(H.NE.ZERO) THEN J=LOG10(ABS(H)) IF(ABS(J).GE.2) THEN ISCALE=ISCALE+J H=H/TEN**J END IF END IF 20 CONTINUE WT=H0*H/DK2/BK1 IF(WT.NE.ZERO) THEN ITEST=LOG10(ABS(WT))+ISCALE IF(ITEST.GE.-NEXP) THEN WT=WT*TEN**ISCALE ELSE WT=ZERO END IF END IF RETURN END IF C General case BK2=ZERO BK1=ZERO RK2=ZERO RK1=ZERO CALL RECUR(N,CK,DK,EK) CALL RECUR(N+1,CK1,DK1,EK1) H=ONE ISCALE=0 DO 10 K=N,1,-1 IF(K.GE.M) THEN RS=T(K-M)/H C Scale and adjust for possible overflow/underflow IF(ISCALE.GT.NEXP) THEN RS=ZERO ELSE RS=RS/TEN**ISCALE END IF ELSE RS=ZERO END IF BB=RS+(DK+XNODE*CK)*BK1+EK1*BK2 BK2=BK1 BK1=BB CALL RECUR(K-1,CKM1,DKM1,EKM1) IF(N.NE.M) H=-H*CK/EK/CKM1 BB=BB*CKM1 WT=BB+(DKM1+XNODE*CKM1)*RK1+EK*RK2 RK2=RK1 RK1=WT CK1=CK DK1=DK EK1=EK CK=CKM1 DK=DKM1 EK=EKM1 IF(BK1.NE.ZERO) THEN J=LOG10(ABS(BK1)) IF(ABS(J).GT.2) THEN C Scale to control overflow/underflow ISCALE=ISCALE+J SCALE=TEN**J BK1=BK1/SCALE BK2=BK2/SCALE RK1=RK1/SCALE RK2=RK2/SCALE END IF END IF 10 CONTINUE WT=H0*BB/WT RETURN END SUBROUTINE TRANSF(T,M,N,RECUR,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION T(0:*) INTEGER M,N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Scales the polynomial expansion: C SUM (M to N) TI*P(I,X). C with respect to the moments HI of the orthogonality weight function C giving the expansion: C H0* SUM (M to N) (TI/HI)*P(I,X). C or C (1/H0)* SUM (M to N) (TI*HI)*P(I,X). C depending on the value of IFLAG. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C C T = Array holding the coefficients TI of the polynomial expansion C to be scaled and expressed as: C SUM (I=M to N) TI*P(I,X) C T(I-M) holds the value of TI. T(*) should be declared to C have at least N-M+1 elements in the calling program. C M = Lower limit to the expansion of T C N = Upper limit to expansion of T C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C IFLAG = 0, if coefficient TI is to be replaced by TI*(H0/HI). C IFLAG = 1, if coefficient TI is to be replaced by TI*(HI/H0). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C T = Array holding the coefficients of the scaled polynomial C expansion. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0D0) H=ONE DO 10 K=0,N CALL RECUR(K,CK,DK,EK) IF(K.NE.0) H=-CKM1/CK*EK*H IF(K.GE.M) THEN IF(IFLAG.EQ.0) THEN T(K-M)=T(K-M)/H ELSE T(K-M)=T(K-M)*H END IF END IF CKM1=CK 10 CONTINUE RETURN END SUBROUTINE BAIR(N,POLIN,POLOUT,A0,A1,A2,RECUR,IDIGIT,ERRVAL,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION POLIN(0:*),POLOUT(0:*),A0,A1,A2,ERRVAL INTEGER N,IDIGIT,IFAIL EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Generalised Bairstow root extraction for polynomial C SUM(I=0 to N) POLIN(I)*P(I,X) C Calculates root as quadratic factor, C A2*P(2,X)-A1*P(1,X)-A0*P(0,X) C where P(I,X) is a general orthogonal polynomial of degree I C C (Reference: Golub & Robertson, Comm.ACM.,10,1967,371-373). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Degree of input polynomial POLIN C POLIN = Coefficients of polynomial of degree N whose quadratic C factor is to be found, i.e. C POLIN = SUM(I=0 to N) POLIN(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C A0,A1,A2 = Estimated quadratic factors C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for full description. C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least signific C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C POLOUT = Coefficients of the deflated polynomial of degree N-2 with C quadratic factor removed, i.e. C POLOUT = SUM(I=0 to N-2) POLOUT(I)*P(I,X) C This array should be declared to have at least N-1 elements C in the calling program. C A0,A1,A2 = Calculated coefficients of the quadratic factor C IFAIL = 0 Quadratic factor found. C = 1 Convergence not achieved after 50 iterations. C ERRVAL = Mean value of the correction to the coefficients of C the quadratic factor. May be used as a measure of the C root accuracy when convergence is not achieved. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,HALF=0.5D0,ONE=1.0D0,TEN=10.0D0) IFAIL=0 ITER=50 ERRVAL=ZERO C Special cases IF(N.EQ.1) THEN A0=-POLIN(0) A1=-POLIN(1) A2=ZERO RETURN END IF IF(N.EQ.2) THEN A0=-POLIN(0) A1=-POLIN(1) A2= POLIN(2) RETURN END IF C Estimated ALPHA & BETA TOL=TEN**(-MAX(1,IDIGIT)) ALPHA=A1/A2 BETA=A0/A2 10 ITER=ITER-1 IF(ITER.LT.0) THEN IFAIL=1 GOTO 20 END IF CALL QFACT(N,POLIN,POLOUT,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) SCALE=MAX(ABS(AB),ABS(BB)) F1=AB/SCALE F2=BB/SCALE DELTA=(B*F1-A*F2)/(AA*F2-BA*F1) SCALE=MAX(ABS(BA),ABS(AA)) F1=BA/SCALE F2=AA/SCALE EPS=(A*F1-B*F2)/(BB*F2-AB*F1) ALPHA=ALPHA+DELTA BETA=BETA+EPS C Test for convergence C Stop if correction is less than 1/TOL times the smallest machine C relative error. IF(ABS(ALPHA)+TOL*ABS(DELTA).NE.ABS(ALPHA) * .OR. ABS(BETA)+TOL*ABS(EPS).NE.ABS(BETA)) GOTO 10 C Final iteration to tidy up result CALL QFACT(N,POLIN,POLOUT,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) SCALE=MAX(ABS(AB),ABS(BB)) F1=AB/SCALE F2=BB/SCALE DELTA=(B*F1-A*F2)/(AA*F2-BA*F1) SCALE=MAX(ABS(BA),ABS(AA)) F1=BA/SCALE F2=AA/SCALE EPS=(A*F1-B*F2)/(BB*F2-AB*F1) ALPHA=ALPHA+DELTA BETA=BETA+EPS 20 A0=BETA A1=ALPHA A2=ONE ERRVAL=HALF*(ABS(EPS)+ABS(DELTA)) RETURN END SUBROUTINE QFACT(N,GAMMA,DELTA,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION GAMMA(0:*),DELTA(0:*),ALPHA,BETA,A,B,AA,AB,BA,BB INTEGER N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Divide the polynomial SUM(I=0 to N) GAMMA(I)*P(I,X) C by the quadratic factor, P(2,X)-ALPHA*P(1,X)-BETA*P(0,X) C giving the quotient SUM(I=0 to N-2) DELTA(I)*P(I,X) C and remainder A*P(1,X)+B*P(0,X) where P(I,X) is the orthogonal C polynomial of degree I defined by RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Degree of GAMMA C GAMMA = Polynomial to be divided by quadratic factor C ALPHA,BETA = Coefficients of quadratic factor C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C DELTA = Quotient polynomial of degreee N-2 C A,B = Remainder coefficients C AA = Partial of A with respect to ALPHA C AB = Partial of A with respect to BETA C BA = Partial of B with respect to ALPHA C BB = Partial of B with respect to BETA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,ONE=1.0D0) C Initialise coefficients DNP2=ZERO DNP1=ZERO DN =ZERO DNM1=ZERO C Partial coefficients wrt ALPHA ADNP2=ZERO ADNP1=ZERO ADN =ZERO ADNM1=ZERO C Partial coefficients wrt BETA BDNP2=ZERO BDNP1=ZERO BDN =ZERO BDNM1=ZERO C C Scaling parameters SN1=ONE SN2=ONE SN3=ONE SN4=ONE CALL RECUR(0,C0,D0,E0) CALL RECUR(1,C1,D1,E1) CALL RECUR(2,C2,D2,E2) CALL RECUR(3,C3,D3,E3) R0=-C0*E1/C1 R1=-C0*E2/C2 R2=-C0*E3/C3 VM1=D0-C0*D1/C1 VM2=D0-C0*D2/C2 W0=-R1*E1 W1=-C1*R2*E2/C2 V1=D1*R1-C1*VM2*E2/C2-C1*R1*D1/C1 K=N-2 CALL RECUR(K+4,CK4,DK4,EK4) CALL RECUR(K+3,CK3,DK3,EK3) CALL RECUR(K+2,CK2,DK2,EK2) CALL RECUR(K+1,CK1,DK1,EK1) VLK4=C0/CK3 VLK3=C0/CK2 VLK2=C0/CK1 RK3=-C0*EK4/CK4 RK2=-C0*EK3/CK3 VMK3=D0-DK3*VLK4 VMK2=D0-DK2*VLK3 C Extract quadratic factor and find partial derivatives DO 100 K=N-2,0,-1 CALL RECUR(K,CK,DK,EK) VLK1=C0/CK RK1=-C0*EK2/CK2 VMK1=D0-DK1*VLK2 SK2=C1*VLK1*VLK2/C0 TK2=VLK2*(D1-C1*DK2/CK2)+C1*VMK1/CK1 UK2=D1*VMK2+E1-C1*VLK3*EK3/CK3-C1*VMK2*DK2/CK2+C1*RK1/CK1 VK2=D1*RK2-C1*VMK3*EK3/CK3-C1*RK2*DK2/CK2 WK2=-C1*RK3*EK3/CK3 CF1=(ALPHA*VLK2-TK2)/SN1 CF2=(BETA+ALPHA*VMK2-UK2)/SN2 CF3=(ALPHA*RK2-VK2)/SN3 CF4=-WK2/SN4 RS=GAMMA(K+2) 40 D=RS+CF1*DNM1+CF2*DN+CF3*DNP1+CF4*DNP2 DELTA(K)=D/SK2 80 DA=VLK2*DNM1/SN1+VMK2*DN/SN2+RK2*DNP1/SN3 * + CF1*ADNM1+CF2*ADN+CF3*ADNP1+CF4*ADNP2 DB=DN/SN2+CF1*BDNM1+CF2*BDN+CF3*BDNP1+CF4*BDNP2 C Recycle old values SN4=SN3 SN3=SN2 SN2=SN1 SN1=SK2 DNP2=DNP1 DNP1=DN DN=DNM1 DNM1=D ADNP2=ADNP1 ADNP1=ADN ADN=ADNM1 ADNM1=DA BDNP2=BDNP1 BDNP1=BDN BDN=BDNM1 BDNM1=DB CK4=CK3 CK3=CK2 CK2=CK1 CK1=CK DK4=DK3 DK3=DK2 DK2=DK1 DK1=DK EK4=EK3 EK3=EK2 EK2=EK1 EK1=EK VLK4=VLK3 VLK3=VLK2 VLK2=VLK1 RK3=RK2 RK2=RK1 VMK3=VMK2 VMK2=VMK1 100 CONTINUE CF1=ALPHA CF2=BETA+ALPHA*VM1-R1 CF3=ALPHA*R1-V1 CF4=-W1 CF5=ALPHA*R0 RS0=GAMMA(0) RS1=GAMMA(1) DNM1=DNM1/SN1 DN=DN/SN2 DNP1=DNP1/SN3 DNP2=DNP2/SN4 ADNM1=ADNM1/SN1 ADN=ADN/SN2 ADNP1=ADNP1/SN3 ADNP2=ADNP2/SN4 BDNM1=BDNM1/SN1 BDN=BDN/SN2 BDNP1=BDNP1/SN3 BDNP2=BDNP2/SN4 C Remainder A=RS1+CF1*DNM1+CF2*DN+CF3*DNP1+CF4*DNP2 B=RS0+BETA*DNM1+CF5*DN-W0*DNP1 C Partials AA=DNM1+VM1*DN+R1*DNP1+CF1*ADNM1+CF2*ADN+CF3*ADNP1+CF4*ADNP2 AB=DN+CF1*BDNM1+CF2*BDN+CF3*BDNP1+CF4*BDNP2 BA=R0*DN+BETA*ADNM1+CF5*ADN-W0*ADNP1 BB=DNM1+BETA*BDNM1+CF5*BDN-W0*BDNP1 RETURN END SUBROUTINE ROOTS(A0,A1,A2,ZREAL1,ZIMAG1,ZREAL2,ZIMAG2,RECUR,INFO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A0,A1,A2,ZREAL1,ZIMAG1,ZREAL2,ZIMAG2 INTEGER INFO EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the roots corresponding to the quadratic factor C A2*P(2,X)-A1*P(1,X)-A0*P(0,X) C where P(I,X) is a general orthogonal polynomial of degree I C defined by the recurrence calculated by subroutine RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C A0,A1,A2 = Coefficients of quadratic factor C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for full description. C C Output parameters: C ZREAL1 = Real part of root 1 C ZIMAG1 = Imaginary part of root 1 C ZREAL2 = Real part of root 2 C ZIMAG2 = Imaginary part of root 2 C INFO = 0 Two roots found C = 1 One root only (A2=0) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0,HALF=0.5D0,FOUR=4.0D0) INFO=0 C CALL RECUR(0,C0,D0,E0) IF(A2.EQ.ZERO) THEN ZREAL1=-(A0+A1*D0)/A1/C0 ZREAL2=ZERO ZIMAG1=ZERO ZIMAG2=ZERO INFO=1 RETURN END IF CALL RECUR(1,C1,D1,E1) AA=-C0*C1*A2 BB=-A2*(C0*D1+D0*C1)+C0*A1 CC=-D0*D1*A2-E1*A2+A0+A1*D0 Z=BB*BB-FOUR*AA*CC ZR=SQRT(ABS(Z)) IF(Z.GE.ZERO) THEN ZIMAG1=ZERO ZIMAG2=ZERO ZREAL1=HALF*(-BB-SIGN(ZR,BB))/AA ZREAL2=CC/AA/ZREAL1 ELSE ZREAL1=-HALF*BB/AA ZREAL2=ZREAL1 ZIMAG1=HALF*ZR/AA ZIMAG2=-ZIMAG1 END IF END SUBROUTINE NEWTON(T,N,XNODE,RECUR,IDIGIT,DELTA,ERRVAL,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION T(0:*),DELTA(0:*) DOUBLE PRECISION XNODE,ERRVAL INTEGER N,IDIGIT EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Applies Newton's method to find a single root of the C polynomial T expressed as: C T = SUM (I=0 to N) T(I)*P(I,X) C where P(I,X) are the orthogonal polymonials whose recurrence C relation is defined by RECUR. C C The value of T is found from the remainder when T is divided C by (X-XNODE). The derivative (of the remainder) is C calculated simultaneously. The deflated polynomial C DELTA = SUM (I=0 to N-1) DELTA(I)*P(I,X) C is also computed. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C T = Polynomial whose roots define the nodes of the quadrature rul C and expressed as: C T = SUM (I=0 to N) T(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C N = Degree of the expansion of T. C XNODE = Approximation to root C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C IDIGIT = Node convergence paramter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C C Output parameters: C C XNODE = Required root. C DELTA = Array holding the coefficients of the deflated polynomial C of degree N-1. This array should be declared to have at C least N elements in the calling program. C ERRVAL = Value of the correction. May be used as a measure of the C root accuracy when convergence is not achieved. C IFAIL = 0, Convergence OK. C = 1, Unsatisfactory convergence after 50 iterations. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (TEN=10.0D0) C ITER=50 TOL=TEN**(-MAX(1,IDIGIT)) 10 ITER=ITER-1 IF(ITER.LT.0) THEN IFAIL=1 ERRVAL=ABS(EPS) RETURN END IF CALL LFACT(T,DELTA,N,XNODE,RECUR,R,DR) EPS=-R/DR XNODE=XNODE+EPS IF(ABS(XNODE)+TOL*ABS(EPS).NE.ABS(XNODE)) GOTO 10 C Final iteration CALL LFACT(T,DELTA,N,XNODE,RECUR,R,DR) EPS=-R/DR XNODE=XNODE+EPS IFAIL=0 ERRVAL=ABS(EPS) 40 RETURN END SUBROUTINE LFACT(GAMMA,DELTA,N,XNODE,RECUR,R,DR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION GAMMA(0:*),DELTA(0:*),XNODE,R,DR INTEGER N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Remove the linear factor (X-XNODE) from the polynomial expansion C SUM(I=0 to N) GAMMA(I) P(I,X) C to give the quotient, C SUM (I=0 to N-1) DELTA(I)*P(I,X). C and the remainder and its derivative at XNODE. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C GAMMA = Polynomial from which factor is to be removed C and expressed as: C GAMMA = SUM (I=0 to N) GAMMA(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C N = Degree of GAMMA. C XNODE = Node to be removed. C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C DELTA = Quotient polynomial expressed as: C DELTA = SUM (I=0 to N-1) DELTA(I)*P(I,X) C This array should be declared to have at least N elements C in the calling program. C R = Remainder from division. C DR = Derivative of R with respect to XNODE. C (-R/DR is the Newton correction). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0D0) C BK1=ZERO BK2=ZERO DBK1=ZERO DBK2=ZERO CALL RECUR(N,CK,DK,EK) CALL RECUR(N+1,CK1,DK1,EK1) DO 10 K=N,0,-1 R=GAMMA(K)+(DK+XNODE*CK)*BK1+EK1*BK2 DR=(DK+XNODE*CK)*DBK1+EK1*DBK2+CK*BK1 BK2=BK1 BK1=R DBK2=DBK1 DBK1=DR IF(K.NE.0) THEN CALL RECUR(K-1,CKM1,DKM1,EKM1) DELTA(K-1)=R*CKM1 END IF EK1=EK CK=CKM1 DK=DKM1 EK=EKM1 10 CONTINUE RETURN END SUBROUTINE EPROD(N,J,COEFF,WORK,LW,RECUR,IFAIL) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION COEFF(*),WORK(LW,2) INTEGER N,J,LW,IFAIL EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the expansion of a product of two orthogonal polynomials C C P(N,X)*P(J,X) = SUM (I=N-J to N+J ) COEFF(I)*P(I,X) C C where J must not exceed N. The orthogonal polynomials are defined C by the recurrence relation calculated by the external C subroutine RECUR. C C For proper initialisation the subroutine must first be called C with J=0 and the required value of N. Subsequent calls must be in C the order J=1,2,,,,,N with the appropriate expansion being C generated from previous values and returned in COEFF(*). The C coefficients of P(N-J,X),...., P(N+J,X) are stored in the array C COEFF(1),...,COEFF(2*J+1). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Highest polynomial degree. Note that after the initial C call with J=0 the value of N in this argument is ignored. C J = Current product of P(J,X) with P(N,X) to be calculated. C Note that the subroutine must be first called with J=0 and C the required largest N. Subsequent calls must be C in the order J=1,2,..,N. C WORK = Matrix work area which must be declared in the calling C program to have dimensions at least (2*J+1) by 2. C The contents of this work area must not be altered between C calls by the calling program. C LW = Leading dimension of WORK in the calling program C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C COEFF = Array holding the calculated coefficients. C This array should be declared to have at least 2*J+1 elemen C in the calling program. C IFAIL = 0 Result OK C = 1 J exceeds N C = 2 J has not been called sequentially C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- INTEGER S,SS,IX(2) SAVE IX,SS,LAST PARAMETER (ZERO=0.0D0,ONE=1.0D0) C IFAIL=0 C Initialise IF(J.EQ.0) THEN IX(1)=1 IX(2)=2 COEFF(1)=ONE WORK(1,2)=ONE LAST=0 SS=N RETURN END IF S=SS C Check that J does not exceed S value IF(S.LT.J) THEN IFAIL=1 RETURN END IF C Check that J is used sequentially IF(LAST.NE.J-1) THEN IFAIL=2 RETURN END IF LAST=J J2=J+J CALL RECUR(J-1,CJ1,DJ1,EJ1) IF(J.EQ.1) THEN DO 20 I=1,J2+1 COEFF(I)=ZERO 20 CONTINUE ELSE DO 25 I=1,J2-3 COEFF(I+2)=WORK(I,IX(1))*EJ1 25 CONTINUE COEFF(1) =ZERO COEFF(2) =ZERO COEFF(J2) =ZERO COEFF(J2+1)=ZERO END IF IBOT=S-J+1 ITOP=S+J-1 DO 30 II=IBOT,ITOP I=II-IBOT+1 CALL RECUR(II,CI,DI,EI) COEFF(I+2)=COEFF(I+2)+(WORK(I,IX(2))/CI)*CJ1 COEFF(I+1)=COEFF(I+1)+WORK(I,IX(2))*(DJ1-(CJ1/CI)*DI) COEFF(I)=COEFF(I)-(WORK(I,IX(2))/CI)*CJ1*EI 30 CONTINUE II=IX(1) IX(1)=IX(2) IX(2)=II DO 35 I=1,J2+1 WORK(I,IX(2))=COEFF(I) 35 CONTINUE RETURN END SUBROUTINE GEFA77(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO DOUBLE PRECISION A(LDA,*) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of LINPACK routine DGEFA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C GEFA77 FACTORS A MATRIX BY GAUSSIAN ELIMINATION. C C ON ENTRY C C A THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C SUBROUTINES AND FUNCTIONS C C BLAS subroutines: DAXPY,DSCAL,IDAMAX C These have been renamed DAXPY7,DSCAL7,IDAMX7 C C INTERNAL VARIABLES C DOUBLE PRECISION T DOUBLE PRECISION ZERO,ONE INTEGER IDAMX7,J,K,KP1,L,NM1 PARAMETER (ZERO=0.0D0,ONE=1.0D0) C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = IDAMX7(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED IF (A(L,K) .EQ. ZERO) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -ONE/A(K,K) CALL DSCAL7(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY7(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. ZERO) INFO = N RETURN END SUBROUTINE GESL77(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(*),JOB DOUBLE PRECISION A(LDA,*),B(*) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of LINPACK routine DGESL C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C GESL77 SOLVES THE SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY GEFA77. C C ON ENTRY C C A THE OUTPUT FROM GEFA77. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM GEFA77. C C B DIMENSION (N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 C OR DGEFA HAS SET INFO .EQ. 0 . C C SUBROUTINES AND FUNCTIONS C C BLAS subroutines: DAXPY,DDOT C These have been renamed: DAXPY7,DDOT7 C INTERNAL VARIABLES C INTEGER K,KB,L,NM1 DOUBLE PRECISION DDOT7,T DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY7(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL DAXPY7(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = DDOT7(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + DDOT7(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE DSCAL7(N,DA,DX,INCX) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DSCAL C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C DOUBLE PRECISION DA,DX(*) DOUBLE PRECISION ZERO,ONE INTEGER I,INCX,M,MP1,N,NINCX PARAMETER (ZERO=0.0D0,ONE=1.0D0) C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END INTEGER FUNCTION IDAMX7(N,DX,INCX) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine IDAMAX C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C DOUBLE PRECISION DX(*),DMAX INTEGER I,INCX,IX,N C IDAMX7 = 0 IF( N .LT. 1 ) RETURN IDAMX7 = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = ABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(ABS(DX(IX)).LE.DMAX) GO TO 5 IDAMX7 = I DMAX = ABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DMAX = ABS(DX(1)) DO 30 I = 2,N IF(ABS(DX(I)).LE.DMAX) GO TO 30 IDAMX7 = I DMAX = ABS(DX(I)) 30 CONTINUE RETURN END SUBROUTINE DAXPY7(N,DA,DX,INCX,DY,INCY) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DAXPY C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C DOUBLE PRECISION DX(*),DY(*),DA INTEGER I,INCX,INCY,IXIY,M,MP1,N DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C IF(N.LE.0)RETURN IF (DA .EQ. ZERO) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DDOT7(N,DX,INCX,DY,INCY) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DDOT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C DOUBLE PRECISION DX(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N DOUBLE PRECISION ZERO,ONE PARAMETER (ZERO=0.0D0,ONE=1.0D0) C DDOT7 = ZERO DTEMP = ZERO IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT7 = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT7 = DTEMP RETURN END C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration driver C (Quadruple precision version of package for VAX/VMS) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C The following templates demonstrate the use of the procedure EXTEND C to generate sequences of extended quadrature rules for various weight C functions given the definitions of the orthogonal polynomial 3-term C recurrence relations. ICASE selects the C required sequence as follows: C ICASE = 1 3-point Gauss-Legendre in [-1,1] C ICASE = 2 2-point Gauss-Lobatto in [-1,1] C ICASE = 3 6-point Radau in [-1,1] C ICASE = 4 2-point Gauss-Laguerre in [0,infinity) C ICASE = 5 3-point Gauss-Hermite in (-infinity,infinity) C ICASE = 6 3-point Gauss-Jacobi in [0,1] C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 code C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- IMPLICIT REAL*16 (A-H,O-Z) PARAMETER (LDA=257,LDB=2*LDA+1,IDIGIT=8,NEXP=4931,NTEST=4) PARAMETER (ONE=1.0Q0,TWO=2.0Q0,THREE=3.0Q0) REAL*16 T(0:LDA),EXT(0:LDA) REAL*16 QR(LDA),QI(LDA) REAL*16 PNODES(LDA),WT(LDA) REAL*16 ERR(LDA),TEST(0:NTEST) REAL*16 WORKA(LDA,LDA),WORKB(LDB,3) INTEGER IWORK(LDA) LOGICAL SYMMET,START EXTERNAL RECURA,RECURB,RECURC,RECURD C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Get data from terminal C C Select demonstration 5 WRITE(*,10) 10 FORMAT(' Case ?:') READ(*,*)ICASE ICASE=MAX(1,ICASE) C Select number of iterative extensions to be performed WRITE(*,20) 20 FORMAT(' No. of rules?:') READ(*,*)NSEQ NSEQ=MAX(1,NSEQ) I=1 50 GOTO(100,200,300,400,500,600),ICASE C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 1 - extension of 3-point Gauss-Legendre rule C Recurrence defined by RECURA 100 IF(I.EQ.1) THEN C Generate 3-point Gauss initially from zero point rule, C i.e. no pre-assigned nodes, symmetry exploited. WRITE(*,*) 'Gauss-Legendre 3-point extension' N=0 M=3 M0=0 T(0)=ONE SYMMET=.TRUE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 120 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 120 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 2 - extension of 2-point Lobatto rule C Recurrence defined by RECURA 200 IF(I.EQ.1) THEN C Add one node, Pre-assign -1.0 and 1.0, symmetry exploited. WRITE(*,*) 'Lobatto 2-point extension' N=2 M=1 PNODES(1)=ONE PNODES(2)=-ONE SYMMET=.TRUE. START=.TRUE. ELSE C Add N-1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N-1 START=.FALSE. END IF H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 220 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 220 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 3 - extension of 6-point Radau rule C Recurrence defined by RECURA 300 IF(I.EQ.1) THEN C Add five nodes. Pre-assign -1.0, no symmetry. WRITE(*,*) 'Radau 6-point extension' N=1 M=5 PNODES(1)=-ONE SYMMET=.FALSE. START=.TRUE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF H0=TWO IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURA,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 320 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURA,TEST(K),IERR) 320 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 4 - extension of 2-point Gauss-Laguerre C Recurrence defined by RECURB 400 IF(I.EQ.1) THEN C Generate 2-point rule initially from zero point rule, C i.e. no pre-assigned nodes, no symmetry. WRITE(*,*) 'Gauss-Laguerre 2-point extension' N=0 M=2 M0=0 T(0)=ONE SYMMET=.FALSE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=ONE IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURB,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 420 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURB,TEST(K),IERR) 420 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 5 - extension of 3-point Gauss-Hermite rule C Recurrence defined by RECURC 500 IF(I.EQ.1) THEN C Generate 3-point rule initially from zero point rule, C i.e. no pre-assigned nodes, symmetry exploited. WRITE(*,*) 'Gauss-Hermite 3-point extension' M=3 N=0 M0=0 T(0)=ONE SYMMET=.TRUE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension C Zero moment integral = sqrt(pi) H0=TWO*SQRT(ATAN(ONE)) IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURC,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 520 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURC,TEST(K),IERR) 520 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Demonstration 6 - extension of 3-point Gauss-Jacobi C for weight sqrt(x) in [0,1] C Recurrence defined by RECURD 600 IF(I.EQ.1) THEN C Generate 3-point rule initially from zero point rule, C i.e. no pre-assigned nodes, no symmetry. WRITE(*,*) 'Gauss-Jacobi 3-point extension' M=3 N=0 M0=0 T(0)=ONE SYMMET=.FALSE. START=.FALSE. ELSE C Add N+1 nodes using the pre-assigned nodes defined by the T polynomial C generated in the previous cycle M=N+1 START=.FALSE. END IF C Calculate extension H0=TWO/THREE IDEG=N+2*M-1 CALL EXTEND(N,M,M0,T,RECURD,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QR,QI,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) C Tests IF(IFLAG.EQ.0.OR.IFLAG.EQ.6) THEN DO 620 K=0,MIN(NTEST,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECURD,TEST(K),IERR) 620 CONTINUE END IF GOTO 2000 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Display results 2000 WRITE(*,*)'Iteration',I WRITE(*,*)'Coefficients of expansion whose roots ', * 'are the new nodes:' WRITE(*,3000)(EXT(J),J,J=0,M) C WRITE(*,*)'New nodes' WRITE(*,3100)(QR(K),QI(K),IWORK(K),ERR(K),K=1,NODES) C WRITE(*,*)'New full extended expansion' WRITE(*,3010)(T(J-M0),J,J=M0,N) C WRITE(*,3200)I,N,IFLAG,NODES IF(IFLAG.NE.0.AND.IFLAG.NE.6) THEN WRITE(*,*)'Terminated prematurely - see IFLAG' GOTO 5 END IF C Print rule (positive nodes only if symmetry present) IF(SYMMET) THEN NUM=(N+1)/2 ELSE NUM=N END IF WRITE(*,3300)(J,PNODES(J),WT(J),J=1,NUM) C Display test results DO 2010 K=0,MIN(NTEST,IDEG/2) WRITE(*,3400)K,TEST(K) 2010 CONTINUE IF(IFLAG.EQ.6) THEN WRITE(*,*)'Rule test is unsatisfactory' GOTO 5 END IF I=I+1 IF(I.LE.NSEQ) GOTO 50 GOTO 5 C 3000 FORMAT(D25.16,'*P(',I3,',X)') 3010 FORMAT(D25.16,'*P(',I3,',X)/HI') 3100 FORMAT(21X,'REAL',16X,'IMAGINARY',1X,'FLAG',7X,'ERR'/, * (2D25.16,I5,D10.1)) 3200 FORMAT(' Complete extended rule: STEP=',I2,' POINTS=',I3, * ' IFLAG=',I1,' NODES ADDED=',I3) 3300 FORMAT(2X,'No.',21X,'NODE',19X,'WEIGHT',/,(I5,2D25.16)) 3400 FORMAT(' TEST(',I2,')=',D25.16) END C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C User defined subroutines - Recurrence relations C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- SUBROUTINE RECURA(K,C,D,E) IMPLICIT REAL*16 (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C This is an example of a user supplied subroutine to define the C orthogonal polynomials. C C CALL RECUR(K,C,D,E) gives the coefficients C,D and E such that, C C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameter: C K = Index for relation C Output parameters: C C,D,E = Parameters in the recurrence relation C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0) C Legendre recurrence for [-1,1] C Covers Gauss, Lobatto and Radau F=FLOAT(K+1) C=FLOAT(2*K+1)/F D=ZERO E=-FLOAT(K)/F RETURN END SUBROUTINE RECURB(K,C,D,E) IMPLICIT REAL*16 (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0Q0) C Laguerre recurrence F=FLOAT(K+1) C=-ONE/F D=FLOAT(2*K+1)/F E=-FLOAT(K)/F RETURN END SUBROUTINE RECURC(K,C,D,E) IMPLICIT REAL*16 (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,ONE=1.0Q0,TWO=2.0Q0) C Hermite recurrence C=ONE D=ZERO E=-FLOAT(K)/TWO RETURN END SUBROUTINE RECURD(K,C,D,E) IMPLICIT REAL*16 (A-H,O-Z) INTEGER K C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: See RECURA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0Q0,TWO=2.0Q0,ONEP5=1.5Q0) C Jacobi polynomials in [0,1] C Weight function is (1-x)**(p-q)*x**(q-1) C C Case for weight sqrt(x) i.e. p=3/2 and q=p P=ONEP5 Q=P FK=FLOAT(K) F2K=FK*TWO B=F2K+P BP1=B+ONE X3=(B-TWO)*(B-ONE)*B A1=(B-ONE)*BP1*X3 A2=-X3*(F2K*(FK+P)+Q*(P-ONE)) A3=X3*BP1*(B-ONE) A4=FK*(FK+Q-ONE)*(FK+P-ONE)*(FK+P-Q)*BP1 C=A3/A1 D=A2/A1 E=-A4/A1 RETURN END C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Main algorithm begins here with structure: C -> ASSIGN C C -> GENER ----> EPROD C C EXTEND --------> SOLVE ----> NEWTON/BAIR C v C v -> RSORT C v C CHECK -> WEIGHT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- SUBROUTINE EXTEND(N,M,M0,T,RECUR,SYMMET,START,PNODES,H0,NEXP, * IDIGIT,WT,NODES,QRNODE,QINODE,ERR,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 T(0:*),EXT(0:*) REAL*16 PNODES(*),WT(*),H0 REAL*16 QRNODE(*),QINODE(*),ERR(*) REAL*16 WORKA(0:LDA-1,0:*),WORKB(0:LDB-1,*) INTEGER M0,N,M,IWORK(*),LDA,LDB,NODES,IFLAG,NEXP,IDIGIT LOGICAL SYMMET,START EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Please address queries or comments to: C C T.N.L. Patterson C Department of Applied Mathematics & Theoretical Physics C The Queen's University of Belfast C Belfast, BT9 1NN C N. Ireland C C Tel: International +44 232 245133 Ext. 3792. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the N+M node quadrature rule composed of N pre-assigned nod C together with M nodes chosen optimally to achieve algebraic degree of C precision of at least N+2*M-1. C C The orthogonal system of polynomials associated with the quadrature C weight is defined generally by the recurrence relation specified in th C user supplied subroutine RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C N = Number of pre-assigned nodes (and upper limit to the expansion C Note that if successful this is reset to N+M on completion C (the appropriate value for iterative use). C M = Number of nodes to be optimally added. C M0 = Lower limit to the expansion of T. This is ignored if START C is .TRUE. Note that if successful this is reset to M C on completion (the appropriate value for iterative use). C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M0 to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. Element T(I-M0) holds the C value of TI. C C Note that T is either, C (1) provided explicitly, C (2) generated automatically from the N pre-assigned nodes C given in PNODES(*) (if START is .TRUE.) C or, C (3) generated from a previous call to the subroutine. C This array should be declared to have at least C max(N-M0+1,M+1) elements in the calling program. C C The service subroutine TRANSF can be used to transform C the expansion to the required input form if desired C with the parameter IFLAG set to 1. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C SYMMET = .FALSE. if no advantage is to be taken of symmetry, if any, C about x=0 in the interval of integration and the C orthogonality weight function. Note that if symmetry in C fact does exist setting this parameter to zero will still C produce correct results - only efficiency is effected. C = .TRUE. if the extended rule computations should C exploit symmetry about x=0 in the interval of C integration and the orthogonality weight function. C This reduces the size of the system of linear equations C determining EXT by a factor of about 2 (see WORKA). If C symmetry does not in fact exist erroneous results will be C produced. C START = .TRUE. then the polynomial T is generated to have C the pre-assigned nodes (PNODES) as its roots. C = .FALSE. then the supplied values of the coefficients C of T are used directly (see above). C PNODES = Array holding the pre-assigned nodes. This array should C be declared to have at least N+M elements in the calling prog C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. C NEXP = Largest negative decimal exponent supported on the C computer. (Positive number - typical value 38). C Weights less than approximately 10**(-NEXP) are set to zero C when the Christoffel-Darboux identity is used (N=M). C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C IWORK = Integer working array which should be declared in the C calling program to have at least max(M,N) elements. C On return IWORK provides information on the convergence C of the nodes. See output parameters. C WORKA = Real working matrix which should be declared in the calling C program to have dimension at least max(M+1,N) C by max(M+1,N+1). If SYMMET=.TRUE. (see above) the C dimension can be reduced to max(M/2+1,N) C by max(M/2+1,N+1). C LDA = Number of elements in the leading dimension of WORKA C declared in the calling program. C WORKB = Real working matrix which should be declared in the calling C program to have dimension at least 2*M+1 by 3. C LDB = Number of elements in the leading dimension of WORKB C declared in the calling program C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C PNODES = Ordered array holding the N+M nodes of the extended C quadrature rule made up from the original pre-assigned C nodes and the new optimally extended nodes. These values can C be used in subsequent iterative use of the subroutine. C WT = Array holding the values of the quadrature weights for C the extended rule associated with the nodes held in PNODES. C This array should be declared to have at least N+M elements C in the calling program. C T = Array holding the coefficients TI of the new orthogonal C expansion whose roots are the nodes of the extended quadratur C (that is, the pre-assigned nodes plus the extended nodes) and C is expressed as: C SUM (I=M to N+M) (TI/HI)*P(I,X) C T(I-M) holds the value of TI. C (For definitions see description of input argument T). C This polynomial can be used as input for further extensions. C C The service subroutine TRANSF can be used to remove the C moment factors from the expansion if desired with the C parameter IFLAG set to 0. C M0 = Lower limit defining the new orthogonal expansion T. C (Set to M). C N = Upper limit defining the new orthogonal expansion T. C (Set to the original value of N+M). C NODES = Number of extended nodes found. Normally equals M but see IFL C QRNODE = Array holding the real parts of the extended nodes (1,..,NODE C This array should be declared to have at least M elements C in the calling program. C QINODE = Array holding the imaginary parts of the extended C nodes (1,..,NODES). (Hopefully these values are zero!). C This array should be declared to have at least M elements C in the calling program. C ERR = Array holding a measure of the relative error in the C nodes. This may be inspected if the convergence C error flag has been raised (IFLAG=3) to decide if the nodes C in question are acceptable. (ERR(*) actually gives the mean C last correction to the quadratic factor in the generalised C Bairstow root finder (see BAIR). This should declared in C the calling program to have at least M elements. C EXT = Array holding the coefficients of the polynomial whose C roots are the extended nodes (QRNODES(*),QINODES(*)) and C expressed as: C EXT = SUM (I=0 to M) EXT(I)*P(I,X) C This array should be declared to have at least M+1 elements C in the calling program. C IWORK = Node convergence flags. Elements 1 to NODES give information C on the convergence of the roots of the polynomial EXT C corresponding to each extended node. C Element I = 0 Convergence of I th root satisfactory C Element I = 1 Convergence of I th root unsatisfactory C IFLAG = 0, No error detected C = 1, The linear system of equations defining the polynomial C whose roots are the extended nodes became singular or C very ill-conditioned. (FATAL). C = 2, The linear system of equations used to generate the C polynomial T when START is .TRUE. became singular C or very ill-conditioned. (FATAL). C = 3, Poor convergence has been detected in the calculation C of the roots of EXT (see above) corresponding to the new C nodes or all nodes have not been found (M not equal C to NODES). See also ERR(*) below. C = 4, Possible imaginary nodes detected. C = 5, Value of N and M incompatible for SYMMET=.TRUE. C Both cannot be odd. (FATAL) C = 6, Test of new quadrature rule has failed. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are included and renamed GEFA77 and GESL7 C These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are renamed DSCAL7, IDAMX7, DAXPY7 and DDOT7. C (Quadruple precision versions used for this subprogram) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Changing the precision. C C This is accomplished as follows: C (1) Amend the TYPE statements, C (2) Select an appropriate value for the NEXP argument to EXTEND. C C NOTE: C (a) All constants used are specified in PARAMETER statements at the st C of each subprogram and, C (b) Generic names are used for all function calls. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- INTEGER S C IFLAG=0 NODES=0 IDEG=N+2*M-1 C Look for incompatible values of N and M IF(SYMMET) THEN C Both N and M cannot be odd IF(MOD(N,2).EQ.1.AND.MOD(M,2).EQ.1) THEN IFLAG=5 RETURN END IF END IF C Generate if required the initial T polynomial corresponding to C prescribed pre-assigned nodes IF(START .AND. N.NE.0) THEN CALL ASSIGN(N,PNODES,IWORK,WORKA,LDA,RECUR,T,IERR) M0=0 IF(IERR.NE.0) THEN IFLAG=2 RETURN END IF END IF NLAST=N C Generate extended expansion coefficients and overwrite T CALL GENER(T,M0,N,M,RECUR,SYMMET,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IERR) IF(IERR.NE.0) THEN IFLAG=1 RETURN END IF C Find extended nodes as roots of EXT(*) CALL SOLVE(EXT,M,SYMMET,RECUR,IDIGIT,QRNODE,QINODE, * NODES,ERR,IWORK,WORKB,LDB,IERR) IF(IERR.NE.0) IFLAG=IERR+2 IF(IFLAG.NE.0) RETURN C Accumulate nodes for extended rule DO 10 I=1,M PNODES(NLAST+I)=QRNODE(I) 10 CONTINUE C Re-order CALL RSORT(PNODES,N,1) C Compute weights (only for positive nodes if symmetric) IF(SYMMET) THEN NUM=(N+1)/2 ELSE NUM=N END IF DO 20 I=1,NUM CALL WEIGHT(T,M0,N,PNODES(I),RECUR,H0,NEXP,WT(I)) IF(SYMMET) THEN WT(N-I+1)=WT(I) END IF 20 CONTINUE C Test the new rule DO 30 K=0,MIN(4,IDEG/2) CALL CHECK(N,PNODES,WT,K,H0,RECUR,TEST,IERR) IF(IERR.EQ.1) THEN IFLAG=6 RETURN END IF 30 CONTINUE RETURN END SUBROUTINE CHECK(N,QNODE,WT,K,H0,RECUR,TEST,IERR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 QNODE(*),WT(*),H0,TEST INTEGER N,K,IERR EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Carry out test of the given quadrature rule by computing the C appropriate integral of, C W(X)*P(K,X)*P(K,X) C over the region associated with the weight function W(X) and the C orthogonal polynomials P(K,X). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Number of nodes in the quadrature rule. C QNODE = Array holding the N nodes. C WT = Array holding the N weights. C K = Index of the orthogonal polynomial whose weighted square C is to be integrated. C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. Note that C P(0,X) is arbitrarily taken to be 1.0. C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C C Output parameters: C TEST = Approximate value of the test integral normalised to C unity. Thus, ABS(TEST-1) gives a measure of the C quality of the calculated rule. C IERR = 0, OK. C = 1, Rule quality unsatisfactory C = 2, Invalid values for input arguments C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,ONE=1.0Q0,TOL=0.0000001Q0) IERR=0 IF(K.LT.0 .OR. N.LT.1 .OR. H0.LE.ZERO) THEN IERR=2 RETURN END IF TEST=ZERO DO 30 I=1,N P1=ONE IF(K.EQ.0) GOTO 20 P0=ZERO X=QNODE(I) C Calculate integrand DO 10 J=0,K-1 CALL RECUR(J,CJ,DJ,EJ) P=(CJ*X+DJ)*P1+EJ*P0 P0=P1 P1=P 10 CONTINUE 20 TEST=TEST+P1*P1*WT(I) 30 CONTINUE TEST=TEST/H0 IF(K.EQ.0) RETURN C Calculate exact value CALL RECUR(0,P,P0,P1) DO 70 J=1,K CALL RECUR(J,CJ,DJ,EJ) P=-P*EJ 70 CONTINUE C Normalise result to unity TEST=TEST*CJ/P C Test for rule quality IF(ABS(TEST-ONE).GT.TOL) IERR=1 RETURN END SUBROUTINE ASSIGN(N,PNODES,IWORK,WORK,LDW,RECUR,T,IERR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 PNODES(*),WORK(0:LDW-1,0:*),T(0:*) INTEGER N,LDW,IERR,IWORK(*) EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Generate the initial polynomial T whose roots are the required C pre-assigned nodes C C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C N = Number of pre-assigned nodes to be used to generate T. C PNODES = Array holding N pre-assigned nodes to be be used to C generate T. C IWORK = Integer working array which should be declared in the C calling program to have at least N elements. C WORK = Real working matrix which should be declared in C the calling program to have dimension at least N by C N+1. C LDW = Number of elements in the leading dimension of WORK C declared in the calling program C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C T = Array holding the coefficients of the polynomial whose C roots define the pre-assigned nodes of the quadrature C rule and expressed as: C H0* SUM (I=0 to N) T(I)/HI*P(I,X) C T(I) holds the value of TI. C This array should be declared to have at least N+1 elements C in the calling program. C IERR = 0, No error detected C = 1, The linear system of equations used to generate the C polynomial T became singular or very ill-conditioned. (FAT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C External library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are used and renamed GEFA77 and GESL77. C Note -- These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are not renamed. C (Quadruple precision versions used for this subprogram) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) IERR=0 C Set up the linear system of equations DO 20 J=1,N X=PNODES(J) P0=ZERO P1=ONE P=P1 DO 10 K=0,N WORK(J-1,K)=P CALL RECUR(K,C0,D0,E0) P=(C0*X+D0)*P1+E0*P0 P0=P1 P1=P 10 CONTINUE 20 CONTINUE C Solve linear system CALL GEFA77(WORK,LDW,N,IWORK,INFO) IF(INFO.NE.0) THEN IERR=1 RETURN END IF CALL GESL77(WORK,LDW,N,IWORK,WORK(0,N),0) DO 30 J=0,N-1 T(J)=-WORK(J,N) 30 CONTINUE T(N)=ONE C Weight with moments CALL TRANSF(T,0,N,RECUR,1) RETURN END SUBROUTINE GENER(T,M0,N,M,RECUR,SYMMET,EXT, * IWORK,WORKA,LDA,WORKB,LDB,IFLAG) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 WORKA(0:LDA-1,0:*),WORKB(0:LDB-1,*) REAL*16 T(0:*),EXT(0:*) INTEGER M0,N,M,IWORK(*),LDA,LDB,IFLAG LOGICAL SYMMET C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Given N pre-assigned quadrature nodes defined as the roots of the C polynomial expansion, C SUM (I=M0 to N) (TI/HI)*P(I,X) C calculate the polynomial expansion, C SUM (I=0 to M) SI*P(I,X) C whose roots are the M optimal nodes and new expansion C SUM (I=M to N+M) (RI/HI)*P(I,X) C whose roots are to the N+M nodes of the full extended quadrature rule. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M0 to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. T(I-M0) holds the C value of TI. This array should be declared to have at least C max(N-M0+1,M+1) elements in the calling program. C M0 = Lower limit to the expansion of T. C N = Upper limit to expansion of T. C M = Number of nodes to be optimally added. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C SYMMET= .FALSE. if no advantage is to be taken of symmetry, if any, C about x=0 in the interval of integration and the C orthogonality weight function. Note that if symmetry in C fact does exist setting this parameter to zero will still C produce correct results - only efficiency is effected. C = .TRUE. if the extended rule computations should C exploit symmetry about x=0 in the interval of C integration and the orthogonality weight function. C This reduces the size of the system of linear equations C determining EXT by a factor of about 2 (see WORKA). If C symmetry does not in fact exist erroneous results will be C produced. C IWORK = Integer working array which should be declared in the C calling program to have at least M elements. C WORKA = Real working matrix which should be declared in the calling C program to have dimension at least M+1 by max(M+1,N+1). C If SYMMET=.TRUE. (see above) the dimension can be reduced to C M/2+1 by max(M/2+1,N/2+1). C LDA = Number of elements in the leading dimension of WORKA C declared in the calling program C WORKB = Real working matrix which should be declared in the calling C program to have at dimension at least 2*M+1 by 3. C LDB = Number of elements in the leading dimension of WORKB C declared in the calling program. C C Output parameters: C T = Array holding the coefficients of the new orthogonal C expansion whose roots are the nodes of the extended quadratur C (that is the pre-assigned nodes plus the extended nodes). C It is expressed as: C SUM (I=M to N+M) (TI/HI)*P(I,X) C where N and M have their original values. T(I-M) holds C the value of TI. See input argument of T for definitions. C M0,N = Lower and upper limits defining the new orthogonal expansion C EXT = Array holding the coefficients of the polynomial whose C roots are the new extended nodes and expressed as: C EXT = SUM (I=0 to M) EXT(I)*P(I,X) C IFLAG = 0, No error detected C = 1, The linear system of equations defining the polynomial C whose roots are the extended nodes became singular or C very ill-conditioned. (FATAL). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C External library routines called: LINPACK - DGEFA, DGESL C FORTRAN-77 versions of these are used and renamed GEFA77 and GESL77. C Note -- These call the BLAS routines DSCAL, IDAMAX, DAXPY and DDOT C which are not renamed. C (Quadruple precision versions used for this subprogram) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- LOGICAL NEVEN,MSODD,MISS EXTERNAL RECUR INTEGER S PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) IFLAG=0 C Look for trivial case IF(N.EQ.0) THEN DO 10 I=0,M-1 EXT(I)=ZERO 10 CONTINUE EXT(M)=ONE T(0)=ONE N=M M0=M RETURN END IF C General case NEVEN=MOD(N,2).EQ.0 NM=N+M C Form matrix 20 DO 60 S=0,M MSODD=MOD(M+S,2).EQ.1 IF(NEVEN.AND.MSODD.AND.SYMMET) GOTO 60 DO 50 J=0,S CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IF(MOD(N+S+J,2).EQ.1.AND.SYMMET) GOTO 50 IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 40 DO 30 I=IBOT,ITOP SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 30 CONTINUE 40 IF(.NOT.SYMMET) THEN WORKA(S,J)=SUM WORKA(J,S)=SUM GOTO 50 END IF IF(NEVEN) THEN WORKA(S/2,J/2)=SUM WORKA(J/2,S/2)=SUM ELSE IF(MSODD) THEN WORKA(S/2,J/2)=SUM ELSE WORKA(J/2,S/2)=SUM END IF END IF 50 CONTINUE 60 CONTINUE NEQ=M IF(SYMMET) NEQ=M/2 C Solve for expansion coefficients CALL GEFA77(WORKA,LDA,NEQ,IWORK,INFO) IF(INFO.NE.0) THEN IFLAG=1 RETURN END IF CALL GESL77(WORKA,LDA,NEQ,IWORK,WORKA(0,NEQ),0) C Store expansion coefficients DO 70 J=0,NEQ-1 EXT(J)=-WORKA(J,NEQ) 70 CONTINUE EXT(NEQ)=ONE C Calculate new T polynomial IF(SYMMET) GOTO 160 C C Non-symmetric case DO 140 S=M,NM IF(S.EQ.M) GOTO 120 DO 110 J=0,M CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 100 DO 90 I=IBOT,ITOP IR=I-IREF SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 90 CONTINUE 100 WORKA(M,J)=SUM 110 CONTINUE 120 SUM=ZERO DO 130 I=0,M SUM=SUM+EXT(I)*WORKA(M,I) 130 CONTINUE WORKA(M-1,S-M)=SUM 140 CONTINUE C Overwrite old values of T DO 150 I=0,N T(I)=WORKA(M-1,I) 150 CONTINUE GOTO 250 C C Symmetric case 160 DO 210 S=M,NM IF(MOD(N+M+S,2).EQ.1) GOTO 210 DO 190 J=0,M CALL EPROD(S,J,WORKB(0,1),WORKB(0,2),LDB,RECUR,IFAIL) IF(MOD(N+S+J,2).EQ.1) GOTO 190 IREF=S-J ITOP=MIN(N,J+S) IBOT=MAX(M0,IREF) SUM=ZERO IF(IBOT.GT.ITOP) GOTO 180 DO 170 I=IBOT,ITOP IR=I-IREF SUM=SUM+T(I-M0)*WORKB(I-IREF,1) 170 CONTINUE 180 WORKA(NEQ,J/2)=SUM 190 CONTINUE SUM=ZERO DO 200 I=0,NEQ SUM=SUM+EXT(I)*WORKA(NEQ,I) 200 CONTINUE WORKA(NEQ-1,(S-M)/2)=SUM 210 CONTINUE C Overwrite old values of T in full unsymmetric form IC=N/2 MISS=.TRUE. DO 220 J=N,0,-1 MISS=.NOT.MISS IF(MISS) THEN T(J)=ZERO ELSE T(J)=WORKA(NEQ-1,IC) IC=IC-1 END IF 220 CONTINUE C Convert EXT to full unsymmetric form WORKB(M,1)=ONE IC=NEQ-1 MISS=.FALSE. DO 230 J=M-1,0,-1 MISS=.NOT.MISS IF(MISS) THEN WORKB(J,1)=ZERO ELSE WORKB(J,1)=EXT(IC) IC=IC-1 END IF 230 CONTINUE DO 240 J=0,M EXT(J)=WORKB(J,1) 240 CONTINUE C Scale new T polynomial 250 PMAX=ZERO DO 260 I=0,N PMAX=MAX(PMAX,ABS(T(I))) 260 CONTINUE DO 270 I=0,N T(I)=T(I)/PMAX 270 CONTINUE N=NM M0=M RETURN END SUBROUTINE SOLVE(EXT,M,SYMMET,RECUR,IDIGIT,QRNODE,QINODE, * NODES,ERR,ICHECK,WORK,LDW,IERR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 EXT(0:*),WORK(0:LDW-1,*),ERR(*) REAL*16 QRNODE(*),QINODE(*) INTEGER M,NODES,LDW,IERR,ICHECK(*),IDIGIT LOGICAL SYMMET,RESET EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculate the roots of the orthogonal polynomial expansion C expressed as, C SUM (I=0 to M) EXT(I)*P(I,X) C where the array EXT holds the appropriate coefficients. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C EXT = Array holding the coefficients of the polynomial whose C roots are required (nodes of the quadrature rule) C and expressed as: C SUM (I=0 to M) EXT(I)*P(I,X) C The recurrence relation for the orthogonal polynomials C P(I,X) is defined by the subroutine RECUR. C This array should be declared to have at least M+1 elements C in the calling program. C M = Upper limit to expansion EXT (polynomial degree). C SYMMET = .FALSE. if no advantage can be taken of symmetry C about x=0 in the interval of integration and the C orthogonality weight function. C = .TRUE. if symmetry exists about x=0 in the interval of C integration and the orthogonality weight function. C RECUR = Name of user supplied subroutine which defines the orthogonal C polynomials. Given K, CALL RECUR(K,C,D,E) gives C the coefficients C,D and E such that, C P(K+1,X)=(C*X+D)*P(K,X)+E*P(K-1,X) C The parameters are defined as follows: C K = Index C C,D,E = Parameters in the recurrence relation C (functions of K) C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C WORK = Real working matrix which should be declared in the calling C program to have dimension at least M+1 by 2. C LDW = Number of elements in the leading dimension of WORK C declared in the calling program C C Output parameters: C C QRNODE = Array holding the real parts of the roots fo EXT (1,..,NODES) C QINODE = Array holding the imaginary parts of the roots of EXT (1,..,N C (hopefully these values are zero!). C NODES = Number of extended nodes found. Normally equals M but see IER C ICHECK = Root convergence flags. Elements 1 to NODES give information C on the convergence of the roots of the polynomial EXT. C Element I = 0 Convergence of I th root satisfactory C Element I = 1 Convergence of I th root unsatisfactory C This array should be declared to have at least M elements C in the calling program. C IERR = 0, No error detected C = 1, Possible imaginary nodes detected. C = 2, Poor convergence has been detected in the calculation C of the roots of EXT (see above) or all roots have not C been found (M not equal to NODES). See also ERR(*) below. C ERR = Array holding a measure of the relative error in the C roots. This may be inspected if the convergence C error flag has been raised (IERR=2) to decide if the roots C in question are acceptable. (ERR(*) actually gives the mean C last correction to the quadratic factor in the generalised C Bairstow root finder (see BAIR). This should declared in C the calling program to have at least M elements. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,HALF=0.5Q0,ONE=1.0Q0, * TWO=2.0Q0,VRT1=0.0000001Q0) NODES=0 IERR=0 C If M is odd find and remove initial real root using NEWTON iteration C Set WORK(*,1) to polynomial to be processed IF(MOD(M,2).EQ.1) THEN ZR1=VRT1 CALL NEWTON(EXT,M,ZR1,RECUR,IDIGIT,WORK(0,1),ERRVAL,IFAIL) NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR1 QINODE(NODES)=ZERO NROOT=M-1 ELSE DO 10 I=0,M WORK(I,1)=EXT(I) 10 CONTINUE NROOT=M END IF IF(NROOT.EQ.0) GOTO 50 C Find remaining root pairs C Calculate seed approximation for quadratic factor CALL RECUR(0,C0,D0,E0) CALL RECUR(1,C1,D1,E1) RT1=VRT1 RT2=ZERO IF(SYMMET) RT2=-RT1 P1A=C0*RT1+D0 P1B=C0*RT2+D0 P2A=(C1*RT1+D1)*P1A+E1 P2B=(C1*RT2+D1)*P1B+E1 DET=C0*(RT1-RT2) SA1=(P2A-P2B)/DET SA0=(P1A*P2B-P1B*P2A)/DET RESET=.TRUE. C Alternate approximation which introduces small complex component RT1=VRT1 RT2=VRT1 SFA1=(C0*D1+D0*C1)/C0+TWO*C1*RT1 SFA0=D0*D1+E1-D0*SFA1-C0*C1*(RT1*RT1+RT2*RT2) C IP1 points to current deflated polynomial IP1=1 IP2=2 20 IF(RESET) THEN A2=ONE A1=SA1 A0=SA0 RESET=.FALSE. END IF CALL BAIR(NROOT,WORK(0,IP1),WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) IF(IFAIL.NE.0) THEN C Try again with complex components introduced A2=ONE A1=SFA1 A0=SFA0 RESET=.TRUE. CALL BAIR(NROOT,WORK(0,IP1),WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) END IF C Apply Bairstow to full expansion to avoid error accumulation CALL BAIR(M,EXT,WORK(0,IP2),A0,A1,A2, * RECUR,IDIGIT,ERRVAL,IFAIL) C Tidy up the quotient polynomial CALL QFACT(NROOT,WORK(0,IP1),WORK(0,IP2),RECUR,A1,A0, * ZR1,ZR1,ZR1,ZR1,ZR1,ZR1) CALL ROOTS(A0,A1,A2,ZR1,ZI1,ZR2,ZI2,RECUR,INFO) C Record node information NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR1 QINODE(NODES)=ZI1 NODES=NODES+1 ICHECK(NODES)=IFAIL ERR(NODES)=ERRVAL QRNODE(NODES)=ZR2 QINODE(NODES)=ZI2 NROOT=NROOT-2 C Make the deflated polynomial current I=IP1 IP1=IP2 IP2=I IF(NROOT.GT.0) THEN C Scale the deflated polynomial PMAX=ZERO DO 30 I=0,NROOT PMAX=MAX(PMAX,ABS(WORK(I,IP1))) 30 CONTINUE DO 40 I=0,NROOT WORK(I,IP1)=WORK(I,IP1)/PMAX 40 CONTINUE GOTO 20 END IF C Calculation complete - Check for difficulties C Look for poor convergence 50 I=0 DO 60 J=1,NODES I=I+ICHECK(J) 60 CONTINUE IF(NODES.NE.M.OR.I.NE.0) THEN IERR=1 RETURN END IF C Look for possible imaginary nodes DO 70 J=1,NODES IF(QINODE(J).NE.ZERO) THEN IERR=2 RETURN END IF 70 CONTINUE RETURN END SUBROUTINE RSORT(A,N,IFLAG) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 A(*) INTEGER N,IFLAG C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C Carries out a simple ripple sort of A(*). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C A = Array holding the numbers to be sorted C N = Number of elements to be sorted C IFLAG = 0 for ascending sort C = 1 for descending sort C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- LOGICAL DONE,ASCEND ASCEND=IFLAG.EQ.0 C Begin scans DO 30 J=N-1,1,-1 DONE=.TRUE. DO 20 K=1,J IF(ASCEND) THEN K1=K K2=K+1 ELSE K1=K+1 K2=K END IF IF(A(K1).GT.A(K2)) THEN C Exchange elements VAL=A(K1) A(K1)=A(K2) A(K2)=VAL DONE=.FALSE. END IF 20 CONTINUE IF(DONE) RETURN 30 CONTINUE RETURN END SUBROUTINE WEIGHT(T,M,N,XNODE,RECUR,H0,NEXP,WT) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 T(0:*),XNODE,H0,WT INTEGER M,N,NEXP EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the quadrature weight associated with the node XNODE in the C rule whose nodes are defined by the roots of polynomial T. C C The weight is calculated by dividing T by (X-XNODE) to give, C C S(X) = T(X)/(X-XNODE) = SUM (0 to N-1) G(I)*P(I,X). C C S(X) is then divided by (X-XNODE) to give remainder R. C C The weight is finally given by H0*G(0)/R. If N=M the C Christoffel-Darboux identity result is used to reduce extreme C cancellation effects at high degree. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C T = Array holding the coefficients TI of the polynomial whose C roots define the N pre-assigned nodes of the quadrature C rule and expressed as: C SUM (I=M to N) (TI/HI)*P(I,X) C where HI is the integral of W(X)*P(I,X)**2 over the C interval for which orthogonality with respect the weight C W(X) is defined (moment integrals) and P(I,X) is the C orthogonal polynomial of degree I. T(I-M) holds the C value of TI. This array should be declared to have at least C N-M+1 elements in the calling program. C M = Lower limit to the expansion of T C N = Upper limit to expansion of T C XNODE = Node whose weight is required C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C H0 = Integral of the orthogonality weight function over the C interval of integration. Zero moment integral. Note that C P(0,X) is arbitrarily taken to be 1.0 C NEXP = Largest negative decimal exponent supported on the C computer. (Positive number - typical value 38). C Weights less than approximately 10**(-NEXP) are set to zero C when the Christoffel-Darboux identity is used (N=M). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C WT = Weight associated with XNODE. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,ONE=1.0Q0,TEN=10.0Q0) C Check for special case IF(M.EQ.N) THEN C Use Christoffel-Darboux result BK1=ZERO BK2=ONE DK1=ZERO DK2=ZERO ISCALE=0 CALL RECUR(0,H,D0,E0) DO 20 K=0,N-1 CALL RECUR(K,CK,DK,EK) IF(K.NE.0) H=-EK*H BB=(CK*XNODE+DK)*BK2+EK*BK1 DD=(CK*XNODE+DK)*DK2+EK*DK1+CK*BK2 BK1=BK2 BK2=BB DK1=DK2 DK2=DD IF(BK2.NE.ZERO) THEN J=LOG10(ABS(BK2)) IF(ABS(J).GT.2) THEN C Scale to control overflow/underflow ISCALE=ISCALE-2*J SCALE=TEN**J BK2=BK2/SCALE BK1=BK1/SCALE DK1=DK1/SCALE DK2=DK2/SCALE END IF END IF IF(H.NE.ZERO) THEN J=LOG10(ABS(H)) IF(ABS(J).GE.2) THEN ISCALE=ISCALE+J H=H/TEN**J END IF END IF 20 CONTINUE WT=H0*H/DK2/BK1 IF(WT.NE.ZERO) THEN ITEST=LOG10(ABS(WT))+ISCALE IF(ITEST.GE.-NEXP) THEN WT=WT*TEN**ISCALE ELSE WT=ZERO END IF END IF RETURN END IF C General case BK2=ZERO BK1=ZERO RK2=ZERO RK1=ZERO CALL RECUR(N,CK,DK,EK) CALL RECUR(N+1,CK1,DK1,EK1) H=ONE ISCALE=0 DO 10 K=N,1,-1 IF(K.GE.M) THEN RS=T(K-M)/H C Scale and adjust for possible overflow/underflow IF(ISCALE.GT.NEXP) THEN RS=ZERO ELSE RS=RS/TEN**ISCALE END IF ELSE RS=ZERO END IF BB=RS+(DK+XNODE*CK)*BK1+EK1*BK2 BK2=BK1 BK1=BB CALL RECUR(K-1,CKM1,DKM1,EKM1) IF(N.NE.M) H=-H*CK/EK/CKM1 BB=BB*CKM1 WT=BB+(DKM1+XNODE*CKM1)*RK1+EK*RK2 RK2=RK1 RK1=WT CK1=CK DK1=DK EK1=EK CK=CKM1 DK=DKM1 EK=EKM1 IF(BK1.NE.ZERO) THEN J=LOG10(ABS(BK1)) IF(ABS(J).GT.2) THEN C Scale to control overflow/underflow ISCALE=ISCALE+J SCALE=TEN**J BK1=BK1/SCALE BK2=BK2/SCALE RK1=RK1/SCALE RK2=RK2/SCALE END IF END IF 10 CONTINUE WT=H0*BB/WT RETURN END SUBROUTINE TRANSF(T,M,N,RECUR,IFLAG) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 T(0:*) INTEGER M,N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Scales the polynomial expansion: C SUM (M to N) TI*P(I,X). C with respect to the moments HI of the orthogonality weight function C giving the expansion: C H0* SUM (M to N) (TI/HI)*P(I,X). C or C (1/H0)* SUM (M to N) (TI*HI)*P(I,X). C depending on the value of IFLAG. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C C T = Array holding the coefficients TI of the polynomial expansion C to be scaled and expressed as: C SUM (I=M to N) TI*P(I,X) C T(I-M) holds the value of TI. T(*) should be declared to C have at least N-M+1 elements in the calling program. C M = Lower limit to the expansion of T C N = Upper limit to expansion of T C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C IFLAG = 0, if coefficient TI is to be replaced by TI*(H0/HI). C IFLAG = 1, if coefficient TI is to be replaced by TI*(HI/H0). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C T = Array holding the coefficients of the scaled polynomial C expansion. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ONE=1.0Q0) H=ONE DO 10 K=0,N CALL RECUR(K,CK,DK,EK) IF(K.NE.0) H=-CKM1/CK*EK*H IF(K.GE.M) THEN IF(IFLAG.EQ.0) THEN T(K-M)=T(K-M)/H ELSE T(K-M)=T(K-M)*H END IF END IF CKM1=CK 10 CONTINUE RETURN END SUBROUTINE BAIR(N,POLIN,POLOUT,A0,A1,A2,RECUR,IDIGIT,ERRVAL,IFAIL) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 POLIN(0:*),POLOUT(0:*),A0,A1,A2,ERRVAL INTEGER N,IDIGIT,IFAIL EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Generalised Bairstow root extraction for polynomial C SUM(I=0 to N) POLIN(I)*P(I,X) C Calculates root as quadratic factor, C A2*P(2,X)-A1*P(1,X)-A0*P(0,X) C where P(I,X) is a general orthogonal polynomial of degree I C C (Reference: Golub & Robertson, Comm.ACM.,10,1967,371-373). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Degree of input polynomial POLIN C POLIN = Coefficients of polynomial of degree N whose quadratic C factor is to be found, i.e. C POLIN = SUM(I=0 to N) POLIN(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C A0,A1,A2 = Estimated quadratic factors C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for full description. C IDIGIT = Node convergence parameter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least signific C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C POLOUT = Coefficients of the deflated polynomial of degree N-2 with C quadratic factor removed, i.e. C POLOUT = SUM(I=0 to N-2) POLOUT(I)*P(I,X) C This array should be declared to have at least N-1 elements C in the calling program. C A0,A1,A2 = Calculated coefficients of the quadratic factor C IFAIL = 0 Quadratic factor found. C = 1 Convergence not achieved after 50 iterations. C ERRVAL = Mean value of the correction to the coefficients of C the quadratic factor. May be used as a measure of the C root accuracy when convergence is not achieved. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,HALF=0.5Q0,ONE=1.0Q0,TEN=10.0Q0) IFAIL=0 ITER=50 ERRVAL=ZERO C Special cases IF(N.EQ.1) THEN A0=-POLIN(0) A1=-POLIN(1) A2=ZERO RETURN END IF IF(N.EQ.2) THEN A0=-POLIN(0) A1=-POLIN(1) A2= POLIN(2) RETURN END IF C Estimated ALPHA & BETA TOL=TEN**(-MAX(1,IDIGIT)) ALPHA=A1/A2 BETA=A0/A2 10 ITER=ITER-1 IF(ITER.LT.0) THEN IFAIL=1 GOTO 20 END IF CALL QFACT(N,POLIN,POLOUT,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) SCALE=MAX(ABS(AB),ABS(BB)) F1=AB/SCALE F2=BB/SCALE DELTA=(B*F1-A*F2)/(AA*F2-BA*F1) SCALE=MAX(ABS(BA),ABS(AA)) F1=BA/SCALE F2=AA/SCALE EPS=(A*F1-B*F2)/(BB*F2-AB*F1) ALPHA=ALPHA+DELTA BETA=BETA+EPS C Test for convergence C Stop if correction is less than 1/TOL times the smallest machine C relative error. IF(ABS(ALPHA)+TOL*ABS(DELTA).NE.ABS(ALPHA) * .OR. ABS(BETA)+TOL*ABS(EPS).NE.ABS(BETA)) GOTO 10 C Final iteration to tidy up result CALL QFACT(N,POLIN,POLOUT,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) SCALE=MAX(ABS(AB),ABS(BB)) F1=AB/SCALE F2=BB/SCALE DELTA=(B*F1-A*F2)/(AA*F2-BA*F1) SCALE=MAX(ABS(BA),ABS(AA)) F1=BA/SCALE F2=AA/SCALE EPS=(A*F1-B*F2)/(BB*F2-AB*F1) ALPHA=ALPHA+DELTA BETA=BETA+EPS 20 A0=BETA A1=ALPHA A2=ONE ERRVAL=HALF*(ABS(EPS)+ABS(DELTA)) RETURN END SUBROUTINE QFACT(N,GAMMA,DELTA,RECUR,ALPHA,BETA,A,B,AA,AB,BA,BB) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 GAMMA(0:*),DELTA(0:*),ALPHA,BETA,A,B,AA,AB,BA,BB INTEGER N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Divide the polynomial SUM(I=0 to N) GAMMA(I)*P(I,X) C by the quadratic factor, P(2,X)-ALPHA*P(1,X)-BETA*P(0,X) C giving the quotient SUM(I=0 to N-2) DELTA(I)*P(I,X) C and remainder A*P(1,X)+B*P(0,X) where P(I,X) is the orthogonal C polynomial of degree I defined by RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Degree of GAMMA C GAMMA = Polynomial to be divided by quadratic factor C ALPHA,BETA = Coefficients of quadratic factor C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C DELTA = Quotient polynomial of degreee N-2 C A,B = Remainder coefficients C AA = Partial of A with respect to ALPHA C AB = Partial of A with respect to BETA C BA = Partial of B with respect to ALPHA C BB = Partial of B with respect to BETA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C Initialise coefficients DNP2=ZERO DNP1=ZERO DN =ZERO DNM1=ZERO C Partial coefficients wrt ALPHA ADNP2=ZERO ADNP1=ZERO ADN =ZERO ADNM1=ZERO C Partial coefficients wrt BETA BDNP2=ZERO BDNP1=ZERO BDN =ZERO BDNM1=ZERO C C Scaling parameters SN1=ONE SN2=ONE SN3=ONE SN4=ONE CALL RECUR(0,C0,D0,E0) CALL RECUR(1,C1,D1,E1) CALL RECUR(2,C2,D2,E2) CALL RECUR(3,C3,D3,E3) R0=-C0*E1/C1 R1=-C0*E2/C2 R2=-C0*E3/C3 VM1=D0-C0*D1/C1 VM2=D0-C0*D2/C2 W0=-R1*E1 W1=-C1*R2*E2/C2 V1=D1*R1-C1*VM2*E2/C2-C1*R1*D1/C1 K=N-2 CALL RECUR(K+4,CK4,DK4,EK4) CALL RECUR(K+3,CK3,DK3,EK3) CALL RECUR(K+2,CK2,DK2,EK2) CALL RECUR(K+1,CK1,DK1,EK1) VLK4=C0/CK3 VLK3=C0/CK2 VLK2=C0/CK1 RK3=-C0*EK4/CK4 RK2=-C0*EK3/CK3 VMK3=D0-DK3*VLK4 VMK2=D0-DK2*VLK3 C Extract quadratic factor and find partial derivatives DO 100 K=N-2,0,-1 CALL RECUR(K,CK,DK,EK) VLK1=C0/CK RK1=-C0*EK2/CK2 VMK1=D0-DK1*VLK2 SK2=C1*VLK1*VLK2/C0 TK2=VLK2*(D1-C1*DK2/CK2)+C1*VMK1/CK1 UK2=D1*VMK2+E1-C1*VLK3*EK3/CK3-C1*VMK2*DK2/CK2+C1*RK1/CK1 VK2=D1*RK2-C1*VMK3*EK3/CK3-C1*RK2*DK2/CK2 WK2=-C1*RK3*EK3/CK3 CF1=(ALPHA*VLK2-TK2)/SN1 CF2=(BETA+ALPHA*VMK2-UK2)/SN2 CF3=(ALPHA*RK2-VK2)/SN3 CF4=-WK2/SN4 RS=GAMMA(K+2) 40 D=RS+CF1*DNM1+CF2*DN+CF3*DNP1+CF4*DNP2 DELTA(K)=D/SK2 80 DA=VLK2*DNM1/SN1+VMK2*DN/SN2+RK2*DNP1/SN3 * + CF1*ADNM1+CF2*ADN+CF3*ADNP1+CF4*ADNP2 DB=DN/SN2+CF1*BDNM1+CF2*BDN+CF3*BDNP1+CF4*BDNP2 C Recycle old values SN4=SN3 SN3=SN2 SN2=SN1 SN1=SK2 DNP2=DNP1 DNP1=DN DN=DNM1 DNM1=D ADNP2=ADNP1 ADNP1=ADN ADN=ADNM1 ADNM1=DA BDNP2=BDNP1 BDNP1=BDN BDN=BDNM1 BDNM1=DB CK4=CK3 CK3=CK2 CK2=CK1 CK1=CK DK4=DK3 DK3=DK2 DK2=DK1 DK1=DK EK4=EK3 EK3=EK2 EK2=EK1 EK1=EK VLK4=VLK3 VLK3=VLK2 VLK2=VLK1 RK3=RK2 RK2=RK1 VMK3=VMK2 VMK2=VMK1 100 CONTINUE CF1=ALPHA CF2=BETA+ALPHA*VM1-R1 CF3=ALPHA*R1-V1 CF4=-W1 CF5=ALPHA*R0 RS0=GAMMA(0) RS1=GAMMA(1) DNM1=DNM1/SN1 DN=DN/SN2 DNP1=DNP1/SN3 DNP2=DNP2/SN4 ADNM1=ADNM1/SN1 ADN=ADN/SN2 ADNP1=ADNP1/SN3 ADNP2=ADNP2/SN4 BDNM1=BDNM1/SN1 BDN=BDN/SN2 BDNP1=BDNP1/SN3 BDNP2=BDNP2/SN4 C Remainder A=RS1+CF1*DNM1+CF2*DN+CF3*DNP1+CF4*DNP2 B=RS0+BETA*DNM1+CF5*DN-W0*DNP1 C Partials AA=DNM1+VM1*DN+R1*DNP1+CF1*ADNM1+CF2*ADN+CF3*ADNP1+CF4*ADNP2 AB=DN+CF1*BDNM1+CF2*BDN+CF3*BDNP1+CF4*BDNP2 BA=R0*DN+BETA*ADNM1+CF5*ADN-W0*ADNP1 BB=DNM1+BETA*BDNM1+CF5*BDN-W0*BDNP1 RETURN END SUBROUTINE ROOTS(A0,A1,A2,ZREAL1,ZIMAG1,ZREAL2,ZIMAG2,RECUR,INFO) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 A0,A1,A2,ZREAL1,ZIMAG1,ZREAL2,ZIMAG2 INTEGER INFO EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the roots corresponding to the quadratic factor C A2*P(2,X)-A1*P(1,X)-A0*P(0,X) C where P(I,X) is a general orthogonal polynomial of degree I C defined by the recurrence calculated by subroutine RECUR. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C A0,A1,A2 = Coefficients of quadratic factor C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for full description. C C Output parameters: C ZREAL1 = Real part of root 1 C ZIMAG1 = Imaginary part of root 1 C ZREAL2 = Real part of root 2 C ZIMAG2 = Imaginary part of root 2 C INFO = 0 Two roots found C = 1 One root only (A2=0) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0,HALF=0.5Q0,FOUR=4.0Q0) INFO=0 C CALL RECUR(0,C0,D0,E0) IF(A2.EQ.ZERO) THEN ZREAL1=-(A0+A1*D0)/A1/C0 ZREAL2=ZERO ZIMAG1=ZERO ZIMAG2=ZERO INFO=1 RETURN END IF CALL RECUR(1,C1,D1,E1) AA=-C0*C1*A2 BB=-A2*(C0*D1+D0*C1)+C0*A1 CC=-D0*D1*A2-E1*A2+A0+A1*D0 Z=BB*BB-FOUR*AA*CC ZR=SQRT(ABS(Z)) IF(Z.GE.ZERO) THEN ZIMAG1=ZERO ZIMAG2=ZERO ZREAL1=HALF*(-BB-SIGN(ZR,BB))/AA ZREAL2=CC/AA/ZREAL1 ELSE ZREAL1=-HALF*BB/AA ZREAL2=ZREAL1 ZIMAG1=HALF*ZR/AA ZIMAG2=-ZIMAG1 END IF END SUBROUTINE NEWTON(T,N,XNODE,RECUR,IDIGIT,DELTA,ERRVAL,IFAIL) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 T(0:*),DELTA(0:*) REAL*16 XNODE,ERRVAL INTEGER N,IDIGIT EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Applies Newton's method to find a single root of the C polynomial T expressed as: C T = SUM (I=0 to N) T(I)*P(I,X) C where P(I,X) are the orthogonal polymonials whose recurrence C relation is defined by RECUR. C C The value of T is found from the remainder when T is divided C by (X-XNODE). The derivative (of the remainder) is C calculated simultaneously. The deflated polynomial C DELTA = SUM (I=0 to N-1) DELTA(I)*P(I,X) C is also computed. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C T = Polynomial whose roots define the nodes of the quadrature rul C and expressed as: C T = SUM (I=0 to N) T(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C N = Degree of the expansion of T. C XNODE = Approximation to root C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C IDIGIT = Node convergence paramter (integer greater than 0). C An attempt is made to calculate the nodes to the maximum C accuracy possible by the machine precision available. C IDIGIT controls the assessment procedure to take account of C round-off errors and specifies the number of least significan C decimal digits that can be ignored (i.e. attributed C to round-off) in the computed relative error. Typical C value is 5. C C Output parameters: C C XNODE = Required root. C DELTA = Array holding the coefficients of the deflated polynomial C of degree N-1. This array should be declared to have at C least N elements in the calling program. C ERRVAL = Value of the correction. May be used as a measure of the C root accuracy when convergence is not achieved. C IFAIL = 0, Convergence OK. C = 1, Unsatisfactory convergence after 50 iterations. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (TEN=10.0Q0) C ITER=50 TOL=TEN**(-MAX(1,IDIGIT)) 10 ITER=ITER-1 IF(ITER.LT.0) THEN IFAIL=1 ERRVAL=ABS(EPS) RETURN END IF CALL LFACT(T,DELTA,N,XNODE,RECUR,R,DR) EPS=-R/DR XNODE=XNODE+EPS IF(ABS(XNODE)+TOL*ABS(EPS).NE.ABS(XNODE)) GOTO 10 C Final iteration CALL LFACT(T,DELTA,N,XNODE,RECUR,R,DR) EPS=-R/DR XNODE=XNODE+EPS IFAIL=0 ERRVAL=ABS(EPS) 40 RETURN END SUBROUTINE LFACT(GAMMA,DELTA,N,XNODE,RECUR,R,DR) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 GAMMA(0:*),DELTA(0:*),XNODE,R,DR INTEGER N EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Remove the linear factor (X-XNODE) from the polynomial expansion C SUM(I=0 to N) GAMMA(I) P(I,X) C to give the quotient, C SUM (I=0 to N-1) DELTA(I)*P(I,X). C and the remainder and its derivative at XNODE. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Input parameters: C GAMMA = Polynomial from which factor is to be removed C and expressed as: C GAMMA = SUM (I=0 to N) GAMMA(I)*P(I,X) C This array should be declared to have at least N+1 elements C in the calling program. C N = Degree of GAMMA. C XNODE = Node to be removed. C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C DELTA = Quotient polynomial expressed as: C DELTA = SUM (I=0 to N-1) DELTA(I)*P(I,X) C This array should be declared to have at least N elements C in the calling program. C R = Remainder from division. C DR = Derivative of R with respect to XNODE. C (-R/DR is the Newton correction). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- PARAMETER (ZERO=0.0Q0) C BK1=ZERO BK2=ZERO DBK1=ZERO DBK2=ZERO CALL RECUR(N,CK,DK,EK) CALL RECUR(N+1,CK1,DK1,EK1) DO 10 K=N,0,-1 R=GAMMA(K)+(DK+XNODE*CK)*BK1+EK1*BK2 DR=(DK+XNODE*CK)*DBK1+EK1*DBK2+CK*BK1 BK2=BK1 BK1=R DBK2=DBK1 DBK1=DR IF(K.NE.0) THEN CALL RECUR(K-1,CKM1,DKM1,EKM1) DELTA(K-1)=R*CKM1 END IF EK1=EK CK=CKM1 DK=DKM1 EK=EKM1 10 CONTINUE RETURN END SUBROUTINE EPROD(N,J,COEFF,WORK,LW,RECUR,IFAIL) IMPLICIT REAL*16 (A-H,O-Z) REAL*16 COEFF(*),WORK(LW,2) INTEGER N,J,LW,IFAIL EXTERNAL RECUR C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C FORTRAN-77 Version 2.2: March 1987 C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Purpose: C C Calculates the expansion of a product of two orthogonal polynomials C C P(N,X)*P(J,X) = SUM (I=N-J to N+J ) COEFF(I)*P(I,X) C C where J must not exceed N. The orthogonal polynomials are defined C by the recurrence relation calculated by the external C subroutine RECUR. C C For proper initialisation the subroutine must first be called C with J=0 and the required value of N. Subsequent calls must be in C the order J=1,2,,,,,N with the appropriate expansion being C generated from previous values and returned in COEFF(*). The C coefficients of P(N-J,X),...., P(N+J,X) are stored in the array C COEFF(1),...,COEFF(2*J+1). C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Unless indicated otherwise the type of each variable is implied C by the default FORTRAN-77 naming convention. C C Input parameters: C C N = Highest polynomial degree. Note that after the initial C call with J=0 the value of N in this argument is ignored. C J = Current product of P(J,X) with P(N,X) to be calculated. C Note that the subroutine must be first called with J=0 and C the required largest N. Subsequent calls must be C in the order J=1,2,..,N. C WORK = Matrix work area which must be declared in the calling C program to have dimensions at least (2*J+1) by 2. C The contents of this work area must not be altered between C calls by the calling program. C LW = Leading dimension of WORK in the calling program C RECUR = Name of the subroutine which defines the orthogonal C polynomials. See EXTEND for a full description. C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- C Output parameters: C C COEFF = Array holding the calculated coefficients. C This array should be declared to have at least 2*J+1 elemen C in the calling program. C IFAIL = 0 Result OK C = 1 J exceeds N C = 2 J has not been called sequentially C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+- INTEGER S,SS,IX(2) SAVE IX,SS,LAST PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C IFAIL=0 C Initialise IF(J.EQ.0) THEN IX(1)=1 IX(2)=2 COEFF(1)=ONE WORK(1,2)=ONE LAST=0 SS=N RETURN END IF S=SS C Check that J does not exceed S value IF(S.LT.J) THEN IFAIL=1 RETURN END IF C Check that J is used sequentially IF(LAST.NE.J-1) THEN IFAIL=2 RETURN END IF LAST=J J2=J+J CALL RECUR(J-1,CJ1,DJ1,EJ1) IF(J.EQ.1) THEN DO 20 I=1,J2+1 COEFF(I)=ZERO 20 CONTINUE ELSE DO 25 I=1,J2-3 COEFF(I+2)=WORK(I,IX(1))*EJ1 25 CONTINUE COEFF(1) =ZERO COEFF(2) =ZERO COEFF(J2) =ZERO COEFF(J2+1)=ZERO END IF IBOT=S-J+1 ITOP=S+J-1 DO 30 II=IBOT,ITOP I=II-IBOT+1 CALL RECUR(II,CI,DI,EI) COEFF(I+2)=COEFF(I+2)+(WORK(I,IX(2))/CI)*CJ1 COEFF(I+1)=COEFF(I+1)+WORK(I,IX(2))*(DJ1-(CJ1/CI)*DI) COEFF(I)=COEFF(I)-(WORK(I,IX(2))/CI)*CJ1*EI 30 CONTINUE II=IX(1) IX(1)=IX(2) IX(2)=II DO 35 I=1,J2+1 WORK(I,IX(2))=COEFF(I) 35 CONTINUE RETURN END SUBROUTINE GEFA77(A,LDA,N,IPVT,INFO) INTEGER LDA,N,IPVT(*),INFO REAL*16 A(LDA,*) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of LINPACK routine DGEFA C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C GEFA77 FACTORS A MATRIX BY GAUSSIAN ELIMINATION. C C ON ENTRY C C A THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C SUBROUTINES AND FUNCTIONS C C BLAS subroutines: DAXPY,DSCAL,IDAMAX C These have been renamed DAXPY7,DSCAL7,IDAMX7 C C INTERNAL VARIABLES C REAL*16 T REAL*16 ZERO,ONE INTEGER IDAMX7,J,K,KP1,L,NM1 PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING C INFO = 0 NM1 = N - 1 IF (NM1 .LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = IDAMX7(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L C C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED IF (A(L,K) .EQ. ZERO) GO TO 40 C C INTERCHANGE IF NECESSARY C IF (L .EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C C COMPUTE MULTIPLIERS C T = -ONE/A(K,K) CALL DSCAL7(N-K,T,A(K+1,K),1) C C ROW ELIMINATION WITH COLUMN INDEXING C DO 30 J = KP1, N T = A(L,J) IF (L .EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY7(N-K,T,A(K+1,K),1,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N) .EQ. ZERO) INFO = N RETURN END SUBROUTINE GESL77(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(*),JOB REAL*16 A(LDA,*),B(*) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of LINPACK routine DGESL C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C GESL77 SOLVES THE SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY GEFA77. C C ON ENTRY C C A THE OUTPUT FROM GEFA77. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A . C C N INTEGER C THE ORDER OF THE MATRIX A . C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM GEFA77. C C B DIMENSION (N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B , C = NONZERO TO SOLVE TRANS(A)*X = B WHERE C TRANS(A) IS THE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X . C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA . IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF DGECO HAS SET RCOND .GT. 0.0 C OR DGEFA HAS SET INFO .EQ. 0 . C C SUBROUTINES AND FUNCTIONS C C BLAS subroutines: DAXPY,DDOT C These have been renamed: DAXPY7,DDOT7 C INTERNAL VARIABLES C INTEGER K,KB,L,NM1 REAL*16 DDOT7,T REAL*16 ZERO,ONE PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C NM1 = N - 1 IF (JOB .NE. 0) GO TO 50 C C JOB = 0 , SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1 .LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L .EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL DAXPY7(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = 1, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL DAXPY7(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE TRANS(A) * X = B C FIRST SOLVE TRANS(U)*Y = B C DO 60 K = 1, N T = DDOT7(K-1,A(1,K),1,B(1),1) B(K) = (B(K) - T)/A(K,K) 60 CONTINUE C C NOW SOLVE TRANS(L)*X = Y C IF (NM1 .LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + DDOT7(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L .EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END SUBROUTINE DSCAL7(N,DA,DX,INCX) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DSCAL C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C REAL*16 DA,DX(*) REAL*16 ZERO,ONE INTEGER I,INCX,M,MP1,N,NINCX PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C IF(N.LE.0)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END INTEGER FUNCTION IDAMX7(N,DX,INCX) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine IDAMAX C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C REAL*16 DX(*),DMAX INTEGER I,INCX,IX,N C IDAMX7 = 0 IF( N .LT. 1 ) RETURN IDAMX7 = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = ABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(ABS(DX(IX)).LE.DMAX) GO TO 5 IDAMX7 = I DMAX = ABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DMAX = ABS(DX(1)) DO 30 I = 2,N IF(ABS(DX(I)).LE.DMAX) GO TO 30 IDAMX7 = I DMAX = ABS(DX(I)) 30 CONTINUE RETURN END SUBROUTINE DAXPY7(N,DA,DX,INCX,DY,INCY) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DAXPY C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C REAL*16 DX(*),DY(*),DA INTEGER I,INCX,INCY,IXIY,M,MP1,N REAL*16 ZERO,ONE PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C IF(N.LE.0)RETURN IF (DA .EQ. ZERO) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END REAL*16 FUNCTION DDOT7(N,DX,INCX,DY,INCY) C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C This is an FORTRAN-77 adaption of BLAS routine DDOT C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+ C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C REAL*16 DX(*),DY(*),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N REAL*16 ZERO,ONE PARAMETER (ZERO=0.0Q0,ONE=1.0Q0) C DDOT7 = ZERO DTEMP = ZERO IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT7 = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT7 = DTEMP RETURN END This disc contains the files: DEXT22.FOR Double precision version of the package. MEXT22.FOR Quadruple precision version for VAX/VMS. SAMPLE.LOG Output of DEXT22.FOR run under VAX/VMS with the input paramanters set to: Case No. of rules 1 3 Gauss-Legendre 2 3 Gauss-Lobatto 3 3 Gauss-Radau 4 2 Gauss-Laguerre 5 2 Gauss-Hermite 6 3 Gauss-Jacobi Log of output from DEXT22 for: ========================================================= Case No. of rules 1 3 Gauss-Legendre 2 3 Gauss-Lobatto 3 3 Gauss-Radau 4 2 Gauss-Laguerre 5 2 Gauss-Hermite 6 3 Gauss-Jacobi ========================================================= RUN DEXT22 Case ?: 1 No. of rules?: 3 Gauss-Legendre 3-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: 0.0000000000000000D+00*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.0000000000000000D+00*P( 2,X) 0.1000000000000000D+01*P( 3,X) New nodes REAL IMAGINARY FLAG ERR 0.0000000000000000D+00 0.0000000000000000D+00 0 0.0D+00 0.7745966692414834D+00 0.0000000000000000D+00 0 0.0D+00 -0.7745966692414834D+00 0.0000000000000000D+00 0 0.0D+00 New full extended expansion 0.1000000000000000D+01*P( 3,X)/HI Complete extended rule: STEP= 1 POINTS= 3 IFLAG=0 NODES ADDED= 3 No. NODE WEIGHT 1 0.7745966692414834D+00 0.5555555555555555D+00 2 0.0000000000000000D+00 0.8888888888888889D+00 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.9999999999999999D+00 TEST( 2)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: 0.1571268237934905D-01*P( 0,X) 0.0000000000000000D+00*P( 1,X) -0.7407407407407407D+00*P( 2,X) 0.0000000000000000D+00*P( 3,X) 0.1000000000000000D+01*P( 4,X) New nodes REAL IMAGINARY FLAG ERR 0.4342437493468026D+00 0.0000000000000000D+00 0 0.0D+00 -0.4342437493468026D+00 0.0000000000000000D+00 0 0.0D+00 0.9604912687080203D+00 0.0000000000000000D+00 0 0.0D+00 -0.9604912687080203D+00 0.0000000000000000D+00 0 0.0D+00 New full extended expansion 0.0000000000000000D+00*P( 4,X)/HI -0.4444444444444444D+00*P( 5,X)/HI 0.0000000000000000D+00*P( 6,X)/HI 0.1000000000000000D+01*P( 7,X)/HI Complete extended rule: STEP= 2 POINTS= 7 IFLAG=0 NODES ADDED= 4 No. NODE WEIGHT 1 0.9604912687080203D+00 0.1046562260264673D+00 2 0.7745966692414834D+00 0.2684880898683335D+00 3 0.4342437493468026D+00 0.4013974147759622D+00 4 0.0000000000000000D+00 0.4509165386584741D+00 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Iteration 3 Coefficients of expansion whose roots are the new nodes: 0.5772349298067298D-03*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.1407910917434097D-02*P( 2,X) 0.0000000000000000D+00*P( 3,X) 0.2790650691523971D+00*P( 4,X) 0.0000000000000000D+00*P( 5,X) -0.1205049839843540D+01*P( 6,X) 0.0000000000000000D+00*P( 7,X) 0.1000000000000000D+01*P( 8,X) New nodes REAL IMAGINARY FLAG ERR 0.2233866864289669D+00 0.0000000000000000D+00 0 0.2D-17 -0.2233866864289669D+00 0.0000000000000000D+00 0 0.2D-17 0.6211029467372264D+00 0.0000000000000000D+00 0 0.1D-17 -0.6211029467372264D+00 0.0000000000000000D+00 0 0.1D-17 0.8884592328722570D+00 0.0000000000000000D+00 0 0.3D-17 -0.8884592328722570D+00 0.0000000000000000D+00 0 0.3D-17 0.9938319632127550D+00 0.0000000000000000D+00 0 0.1D-17 -0.9938319632127550D+00 0.0000000000000000D+00 0 0.1D-17 New full extended expansion 0.0000000000000000D+00*P( 8,X)/HI -0.2896464423067150D-01*P( 9,X)/HI 0.0000000000000000D+00*P( 10,X)/HI 0.3461310330856246D+00*P( 11,X)/HI 0.0000000000000000D+00*P( 12,X)/HI -0.1000000000000000D+01*P( 13,X)/HI 0.0000000000000000D+00*P( 14,X)/HI 0.7705426896934803D+00*P( 15,X)/HI Complete extended rule: STEP= 3 POINTS= 15 IFLAG=0 NODES ADDED= 8 No. NODE WEIGHT 1 0.9938319632127550D+00 0.1700171962994028D-01 2 0.9604912687080203D+00 0.5160328299707982D-01 3 0.8884592328722570D+00 0.9292719531512452D-01 4 0.7745966692414834D+00 0.1344152552437843D+00 5 0.6211029467372264D+00 0.1715119091363914D+00 6 0.4342437493468026D+00 0.2006285293769890D+00 7 0.2233866864289669D+00 0.2191568584015875D+00 8 0.0000000000000000D+00 0.2255104997982067D+00 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Case ?: 2 No. of rules?: 3 Lobatto 2-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: 0.0000000000000000D+00*P( 0,X) 0.1000000000000000D+01*P( 1,X) New nodes REAL IMAGINARY FLAG ERR 0.0000000000000000D+00 0.0000000000000000D+00 0 0.0D+00 New full extended expansion -0.1000000000000000D+01*P( 1,X)/HI 0.0000000000000000D+00*P( 2,X)/HI 0.4285714285714286D+00*P( 3,X)/HI Complete extended rule: STEP= 1 POINTS= 3 IFLAG=0 NODES ADDED= 1 No. NODE WEIGHT 1 0.1000000000000000D+01 0.3333333333333333D+00 2 0.0000000000000000D+00 0.1333333333333333D+01 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: -0.1428571428571429D+00*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.1000000000000000D+01*P( 2,X) New nodes REAL IMAGINARY FLAG ERR 0.6546536707079772D+00 0.0000000000000000D+00 0 0.0D+00 -0.6546536707079771D+00 0.0000000000000000D+00 0 0.0D+00 New full extended expansion 0.0000000000000000D+00*P( 2,X)/HI -0.1000000000000000D+01*P( 3,X)/HI 0.0000000000000000D+00*P( 4,X)/HI 0.6363636363636364D+00*P( 5,X)/HI Complete extended rule: STEP= 2 POINTS= 5 IFLAG=0 NODES ADDED= 2 No. NODE WEIGHT 1 0.1000000000000000D+01 0.1000000000000000D+00 2 0.6546536707079772D+00 0.5444444444444445D+00 3 0.0000000000000000D+00 0.7111111111111111D+00 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 Iteration 3 Coefficients of expansion whose roots are the new nodes: -0.4746768383132019D-01*P( 0,X) 0.0000000000000000D+00*P( 1,X) -0.1515151515151515D+00*P( 2,X) 0.0000000000000000D+00*P( 3,X) 0.1000000000000000D+01*P( 4,X) New nodes REAL IMAGINARY FLAG ERR 0.3409822659109930D+00 0.0000000000000000D+00 0 0.0D+00 -0.3409822659109930D+00 0.0000000000000000D+00 0 0.0D+00 0.8904055275126688D+00 0.0000000000000000D+00 0 0.2D-17 -0.8904055275126688D+00 0.0000000000000000D+00 0 0.2D-17 New full extended expansion 0.0000000000000000D+00*P( 4,X)/HI -0.3813459268004723D+00*P( 5,X)/HI 0.0000000000000000D+00*P( 6,X)/HI -0.9870129870129869D+00*P( 7,X)/HI 0.0000000000000000D+00*P( 8,X)/HI 0.1000000000000000D+01*P( 9,X)/HI Complete extended rule: STEP= 3 POINTS= 9 IFLAG=0 NODES ADDED= 4 No. NODE WEIGHT 1 0.1000000000000000D+01 0.3064373897707232D-01 2 0.8904055275126688D+00 0.1792626995532074D+00 3 0.6546536707079772D+00 0.2839787780481211D+00 4 0.3409822659109930D+00 0.3342337398164177D+00 5 0.0000000000000000D+00 0.3437620872103631D+00 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Case ?: 3 No. of rules?: 3 Radau 6-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: -0.9090909090909092D-01*P( 0,X) 0.2727272727272728D+00*P( 1,X) -0.4545454545454546D+00*P( 2,X) 0.6363636363636364D+00*P( 3,X) -0.8181818181818182D+00*P( 4,X) 0.1000000000000000D+01*P( 5,X) New nodes REAL IMAGINARY FLAG ERR 0.1240503795052277D+00 0.0000000000000000D+00 0 0.4D-17 0.6039731642527837D+00 0.0000000000000000D+00 0 0.6D-17 -0.3909285467072722D+00 0.0000000000000000D+00 0 0.6D-17 0.9203802858970625D+00 0.0000000000000000D+00 0 0.3D-17 -0.8029298284023471D+00 0.0000000000000000D+00 0 0.3D-17 New full extended expansion 0.1000000000000000D+01*P( 5,X)/HI 0.8461538461538462D+00*P( 6,X)/HI Complete extended rule: STEP= 1 POINTS= 6 IFLAG=0 NODES ADDED= 5 No. NODE WEIGHT 1 0.9203802858970625D+00 0.2015883852534808D+00 2 0.6039731642527837D+00 0.4169013343119077D+00 3 0.1240503795052277D+00 0.5209267831895750D+00 4 -0.3909285467072722D+00 0.4853871884689699D+00 5 -0.8029298284023471D+00 0.3196407532205109D+00 6 -0.1000000000000000D+01 0.5555555555555556D-01 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: 0.1238920331102971D-01*P( 0,X) -0.1115360783931244D-01*P( 1,X) 0.2930949908491627D-01*P( 2,X) -0.2573880815371107D-01*P( 3,X) -0.8051048473617479D+00*P( 4,X) -0.6245367209597196D+00*P( 5,X) 0.7380888520433050D+00*P( 6,X) 0.1000000000000000D+01*P( 7,X) New nodes REAL IMAGINARY FLAG ERR -0.1366416372713493D+00 0.0000000000000000D+00 0 0.9D-18 -0.6212165636007762D+00 0.0000000000000000D+00 0 0.4D-17 0.3761018966748717D+00 0.0000000000000000D+00 0 0.4D-17 -0.8487727816552522D+00 0.0000000000000000D+00 0 0.4D-17 0.7905296024586544D+00 0.0000000000000000D+00 0 0.4D-17 0.9865012767506534D+00 0.0000000000000000D+00 0 0.2D-17 -0.9439342521493506D+00 0.0000000000000000D+00 0 0.2D-17 New full extended expansion -0.8815438503554517D-01*P( 7,X)/HI -0.8472981649987972D-01*P( 8,X)/HI -0.3768076078707027D+00*P( 9,X)/HI -0.4473580103536859D+00*P( 10,X)/HI 0.3851866566831191D+00*P( 11,X)/HI 0.1000000000000000D+01*P( 12,X)/HI 0.5106460045718307D+00*P( 13,X)/HI Complete extended rule: STEP= 2 POINTS= 13 IFLAG=0 NODES ADDED= 7 No. NODE WEIGHT 1 0.9865012767506534D+00 0.3610633344712083D-01 2 0.9203802858970625D+00 0.9797463317823482D-01 3 0.7905296024586544D+00 0.1603549149394082D+00 4 0.6039731642527837D+00 0.2100921244820766D+00 5 0.3761018966748717D+00 0.2427304474580505D+00 6 0.1240503795052277D+00 0.2588076132975782D+00 7 -0.1366416372713493D+00 0.2601366692626690D+00 8 -0.3909285467072722D+00 0.2455064597946105D+00 9 -0.6212165636007762D+00 0.2115512241336131D+00 10 -0.8029298284023471D+00 0.1329820397038908D+00 11 -0.8487727816552522D+00 0.3616137156650629D-01 12 -0.9439342521493506D+00 0.9202712852123144D-01 13 -0.1000000000000000D+01 0.1556904021500967D-01 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Iteration 3 Coefficients of expansion whose roots are the new nodes: -0.2568241533603606D+00*P( 0,X) 0.6273882318514978D+00*P( 1,X) -0.6950114612492999D+00*P( 2,X) 0.3514310704836285D+00*P( 3,X) 0.2414864785768811D+00*P( 4,X) -0.9310510552284522D+00*P( 5,X) 0.1469801497278645D+01*P( 6,X) -0.5792538945795940D+01*P( 7,X) 0.4360460897303821D+01*P( 8,X) -0.1885604119806995D+02*P( 9,X) 0.1887338563566692D+01*P( 10,X) 0.6992128566633852D+02*P( 11,X) -0.1072921826820173D+02*P( 12,X) -0.4426841333562902D+02*P( 13,X) 0.1000000000000000D+01*P( 14,X) New nodes REAL IMAGINARY FLAG ERR 0.9979162551518678D+00 0.0000000000000000D+00 0 0.1D-16 -0.6134267604635425D-02 0.0000000000000000D+00 0 0.1D-16 0.9613067296650284D+00 0.0000000000000000D+00 0 0.7D-17 -0.2655914229442032D+00 0.0000000000000000D+00 0 0.7D-17 0.8632649380482529D+00 0.0000000000000000D+00 0 0.8D-17 0.7034887989758216D+00 0.0000000000000000D+00 0 0.8D-17 0.4941080138619341D+00 0.0000000000000000D+00 0 0.4D-16 0.2520863057103871D+00 0.0000000000000000D+00 0 0.4D-16 -0.7205766157140713D+00 0.0000000000000000D+00 0 0.2D-16 -0.5103354152259295D+00 0.0000000000000000D+00 0 0.2D-16 -0.9812660421998112D+00 0.0000000000000000D+00 0 0.1D-15 -0.8923964308975302D+00 0.0000000000000000D+00 0 0.1D-15 0.2307179086994016D+02 0.0000000000000000D+00 0 0.6D-16 -0.1013669616811477D+01 0.0000000000000000D+00 0 0.6D-16 New full extended expansion 0.5454360425671037D-02*P( 14,X)/HI 0.1275790889143820D-02*P( 15,X)/HI 0.1362458441982789D-01*P( 16,X)/HI 0.1061409558174235D-01*P( 17,X)/HI -0.7800045804654459D-01*P( 18,X)/HI -0.6147544605704888D-01*P( 19,X)/HI -0.2027665031851448D+00*P( 20,X)/HI -0.2559758659655763D+00*P( 21,X)/HI 0.5447045628625731D+00*P( 22,X)/HI 0.1000000000000000D+01*P( 23,X)/HI 0.5759468453196904D-02*P( 24,X)/HI -0.7380381496936674D+00*P( 25,X)/HI -0.3289096723205376D+00*P( 26,X)/HI 0.7335603893034949D-02*P( 27,X)/HI Complete extended rule: STEP= 3 POINTS= 27 IFLAG=6 NODES ADDED= 14 No. NODE WEIGHT 1 0.2307179086994016D+02 -0.1710833635652233D-16 2 0.9979162551518678D+00 0.5761238214206027D-02 3 0.9865012767506534D+00 0.1781397602627082D-01 4 0.9613067296650284D+00 0.3288145280959046D-01 5 0.9203802858970625D+00 0.4904105110221086D-01 6 0.8632649380482529D+00 0.6509033800352046D-01 7 0.7905296024586544D+00 0.8015752890516894D-01 8 0.7034887989758216D+00 0.9361527599467614D-01 9 0.6039731642527837D+00 0.1050581380411569D+00 10 0.4941080138619341D+00 0.1143012662323623D+00 11 0.3761018966748717D+00 0.1213546454148586D+00 12 0.2520863057103871D+00 0.1263455706105022D+00 13 0.1240503795052277D+00 0.1294155830050473D+00 14 -0.6134267604635425D-02 0.1306510180165526D+00 15 -0.1366416372713493D+00 0.1300525005145075D+00 16 -0.2655914229442032D+00 0.1275061882372417D+00 17 -0.3909285467072722D+00 0.1227817936161397D+00 18 -0.5103354152259295D+00 0.1155996181401652D+00 19 -0.6212165636007762D+00 0.1056775942146705D+00 20 -0.7205766157140713D+00 0.9232499519844784D-01 21 -0.8029298284023471D+00 0.6879140063252571D-01 22 -0.8487727816552522D+00 0.3092904062992692D-01 23 -0.8923964308975302D+00 0.5477880738383480D-01 24 -0.9439342521493506D+00 0.4546365767396098D-01 25 -0.9812660421998112D+00 0.2876294640770191D-01 26 -0.1000000000000000D+01 0.5992084797070105D-02 27 -0.1013669616811477D+01 -0.1477098223165217D-03 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.9999999999999862D+00 TEST( 2)= 0.9999999999727656D+00 TEST( 3)= 0.9999999436798130D+00 TEST( 4)= 0.9998820691687861D+00 Rule test is unsatisfactory Case ?: 4 No. of rules?: 2 Gauss-Laguerre 2-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: 0.0000000000000000D+00*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.1000000000000000D+01*P( 2,X) New nodes REAL IMAGINARY FLAG ERR 0.3414213562373095D+01 0.0000000000000000D+00 0 0.0D+00 0.5857864376269049D+00 0.0000000000000000D+00 0 0.0D+00 New full extended expansion 0.1000000000000000D+01*P( 2,X)/HI Complete extended rule: STEP= 1 POINTS= 2 IFLAG=0 NODES ADDED= 2 No. NODE WEIGHT 1 0.3414213562373095D+01 0.1464466094067262D+00 2 0.5857864376269049D+00 0.8535533905932737D+00 TEST( 0)= 0.9999999999999999D+00 TEST( 1)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: 0.6000000000000000D+01*P( 0,X) -0.1500000000000000D+01*P( 1,X) 0.0000000000000000D+00*P( 2,X) 0.1000000000000000D+01*P( 3,X) New nodes REAL IMAGINARY FLAG ERR 0.8396196974011156D+01 0.0000000000000000D+00 0 0.7D-16 0.3019015129944220D+00 -0.1959389276469933D+01 0 0.2D-16 0.3019015129944220D+00 0.1959389276469933D+01 0 0.2D-16 New full extended expansion 0.8125000000000000D+00*P( 3,X)/HI -0.1000000000000000D+01*P( 4,X)/HI 0.4166666666666667D+00*P( 5,X)/HI Complete extended rule: STEP= 2 POINTS= 5 IFLAG=4 NODES ADDED= 3 Terminated prematurely - see IFLAG Case ?: 5 No. of rules?: 2 Gauss-Hermite 3-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: 0.0000000000000000D+00*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.0000000000000000D+00*P( 2,X) 0.1000000000000000D+01*P( 3,X) New nodes REAL IMAGINARY FLAG ERR 0.0000000000000000D+00 0.0000000000000000D+00 0 0.0D+00 0.1224744871391589D+01 0.0000000000000000D+00 0 0.0D+00 -0.1224744871391589D+01 0.0000000000000000D+00 0 0.0D+00 New full extended expansion 0.1000000000000000D+01*P( 3,X)/HI Complete extended rule: STEP= 1 POINTS= 3 IFLAG=0 NODES ADDED= 3 No. NODE WEIGHT 1 0.1224744871391589D+01 0.2954089751509193D+00 2 0.0000000000000000D+00 0.1181635900603677D+01 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: -0.3000000000000000D+01*P( 0,X) 0.0000000000000000D+00*P( 1,X) -0.2000000000000000D+01*P( 2,X) 0.0000000000000000D+00*P( 3,X) 0.1000000000000000D+01*P( 4,X) New nodes REAL IMAGINARY FLAG ERR 0.0000000000000000D+00 -0.4884800789447105D+00 0 0.0D+00 0.0000000000000000D+00 0.4884800789447105D+00 0 0.0D+00 0.2288801605103822D+01 0.0000000000000000D+00 0 0.2D-16 -0.2288801605103822D+01 0.0000000000000000D+00 0 0.2D-16 New full extended expansion 0.0000000000000000D+00*P( 4,X)/HI 0.3809523809523810D+00*P( 5,X)/HI 0.0000000000000000D+00*P( 6,X)/HI 0.1000000000000000D+01*P( 7,X)/HI Complete extended rule: STEP= 2 POINTS= 7 IFLAG=4 NODES ADDED= 4 Terminated prematurely - see IFLAG Case ?: 6 No. of rules?: 3 Gauss-Jacobi 3-point extension Iteration 1 Coefficients of expansion whose roots are the new nodes: 0.0000000000000000D+00*P( 0,X) 0.0000000000000000D+00*P( 1,X) 0.0000000000000000D+00*P( 2,X) 0.1000000000000000D+01*P( 3,X) New nodes REAL IMAGINARY FLAG ERR 0.1647102868965424D+00 0.0000000000000000D+00 0 0.3D-17 0.9008058292716294D+00 0.0000000000000000D+00 0 0.4D-17 0.5498684992164436D+00 0.0000000000000000D+00 0 0.4D-17 New full extended expansion 0.1000000000000000D+01*P( 3,X)/HI Complete extended rule: STEP= 1 POINTS= 3 IFLAG=0 NODES ADDED= 3 No. NODE WEIGHT 1 0.9008058292716294D+00 0.2332816246559149D+00 2 0.5498684992164436D+00 0.3076023676819127D+00 3 0.1647102868965424D+00 0.1257826743288388D+00 TEST( 0)= 0.9999999999999998D+00 TEST( 1)= 0.9999999999999998D+00 TEST( 2)= 0.1000000000000000D+01 Iteration 2 Coefficients of expansion whose roots are the new nodes: 0.8302544967323833D-04*P( 0,X) 0.1808212463080860D-03*P( 1,X) -0.6293935530868694D-01*P( 2,X) 0.0000000000000000D+00*P( 3,X) 0.1000000000000000D+01*P( 4,X) New nodes REAL IMAGINARY FLAG ERR 0.3434982475781609D+00 0.0000000000000000D+00 0 0.2D-17 0.4317458752763439D-01 0.0000000000000000D+00 0 0.2D-17 0.9829837529243084D+00 0.0000000000000000D+00 0 0.8D-17 0.7479904707934259D+00 0.0000000000000000D+00 0 0.8D-17 New full extended expansion -0.1000000000000000D+01*P( 4,X)/HI -0.7540286168214371D+00*P( 5,X)/HI -0.9004851826770350D+00*P( 6,X)/HI 0.5021708298937920D+00*P( 7,X)/HI Complete extended rule: STEP= 2 POINTS= 7 IFLAG=0 NODES ADDED= 4 No. NODE WEIGHT 1 0.9829837529243084D+00 0.4509009780887507D-01 2 0.9008058292716294D+00 0.1136674496440606D+00 3 0.7479904707934259D+00 0.1567638583107757D+00 4 0.5498684992164436D+00 0.1546815936593538D+00 5 0.3434982475781609D+00 0.1161056548660293D+00 6 0.1647102868965424D+00 0.6270539664794700D-01 7 0.4317458752763439D-01 0.1765261572962525D-01 TEST( 0)= 0.1000000000000000D+01 TEST( 1)= 0.1000000000000000D+01 TEST( 2)= 0.1000000000000000D+01 TEST( 3)= 0.1000000000000000D+01 TEST( 4)= 0.1000000000000000D+01 Iteration 3 Coefficients of expansion whose roots are the new nodes: 0.1246085068908996D-07*P( 0,X) 0.5106502089882409D-07*P( 1,X) 0.1857830431079363D-06*P( 2,X) 0.3154953497851658D-04*P( 3,X) 0.3729410486893128D-03*P( 4,X) 0.6566133154308339D-02*P( 5,X) -0.6845782043616039D-01*P( 6,X) -0.1121555778222753D+00*P( 7,X) 0.1000000000000000D+01*P( 8,X) New nodes REAL IMAGINARY FLAG ERR 0.9522115224106390D-01 0.0000000000000000D+00 0 0.2D-16 0.1092616002586163D-01 0.0000000000000000D+00 0 0.2D-16 0.4453534840434267D+00 0.0000000000000000D+00 0 0.5D-17 0.2486386449037345D+00 0.0000000000000000D+00 0 0.5D-17 0.8321177481736263D+00 0.0000000000000000D+00 0 0.7D-18 0.6523614224614975D+00 0.0000000000000000D+00 0 0.7D-18 0.9973759430379304D+00 0.0000000000000000D+00 0 0.1D-16 0.9513731441472555D+00 0.0000000000000000D+00 0 0.1D-16 New full extended expansion -0.1000000000000000D+01*P( 8,X)/HI -0.7219216410522186D+00*P( 9,X)/HI -0.4137098497966033D+00*P( 10,X)/HI 0.6243456378913302D+00*P( 11,X)/HI 0.2379858019599160D+00*P( 12,X)/HI 0.3071737360122939D+00*P( 13,X)/HI -0.3226631546888718D+00*P( 14,X)/HI 0.5874911786061134D-01*P( 15,X)/HI Complete extended rule: STEP= 3 POINTS= 15 IFLAG=0 NODES ADDED= 8 No. NODE WEIGHT 1 0.9973759430379304D+00 0.7250785857083849D-02 2 0.9829837529243084D+00 0.2225370714572648D-01 3 0.9513731441472555D+00 0.4003272827656301D-01 4 0.9008058292716294D+00 0.5689076367068382D-01 5 0.8321177481736263D+00 0.7021992138795099D-01 6 0.7479904707934259D+00 0.7836602038520848D-01 7 0.6523614224614975D+00 0.8066055699011988D-01 8 0.5498684992164436D+00 0.7734675361301653D-01 9 0.4453534840434267D+00 0.6936411431773201D-01 10 0.3434982475781609D+00 0.5805035560699054D-01 11 0.2486386449037345D+00 0.4488514434457869D-01 12 0.1647102868965424D+00 0.3135360123263908D-01 13 0.9522115224106390D-01 0.1889127817157171D-01 14 0.4317458752763439D-01 0.8826110122443174D-02 15 0.1092616002586163D-01 0.2274825544359364D-02 TEST( 0)= 0.1000000000000001D+01 TEST( 1)= 0.1000000000000004D+01 TEST( 2)= 0.1000000000000006D+01 TEST( 3)= 0.1000000000000009D+01 TEST( 4)= 0.1000000000000013D+01