C ALGORITHM 604, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 3, C SEP., 1983, P. 381-383. SUBROUTINE EXTREM(X, N, XINT, NINT, K, IPRINT, MAXIT, EPS2, EPS3, EXT 10 * EPS4, XTILDE, ALPHA, PROD, PZERO, IERR) C THIS SUBROUTINE FINDS A POLYNOMIAL WHICH IS AN EXTREMAL OF THE LINEAR C FUNCTIONAL WHICH IS THE POINT EVALUATION AT THE POINT ZERO. THE C POLYNOMIAL IS CHOSEN FROM THE SPACE OF POLYNOMIALS OF DEGREE LESS C THAN OR EQUAL TO N, WHERE N IS A PARAMETER SET BY THE USER. THE NORM C ON THIS SPACE IS DEFINED TO BE: C C NORM(P) = MAX(ABS(P(T)) : T IS IN S) C C WHERE S IS THE UNION OF THE CLOSED INTERVALS (XINT(1,I),XINT(2,I)), C I = 1, 2, ..., NINT. C C THE REMES ALGORITHM (OUTLINED BELOW) IS USED TO FIND THIS EXTREMAL C POLYNOMIAL. C C C PARAMETERS C C X A REAL ARRAY OF DIMENSION N+1. ON INPUT, X CONTAINS A C STRICTLY INCREASING SEQUENCE OF POINTS IN S WITH X(K)=B AND C X(K+1)=C, WHERE B AND C ARE DEFINED: C B = MAX (T : T IS IN S, AND T .LE. 0) C C = MIN (T : T IS IN S, AND T .GE. 0). C ON RETURN, X CONTAINS THE FINAL SEQUENCE WHICH DEFINES THE C EXTREMAL POLYNOMIAL WHICH IS THE UNIQUE POLYNOMIAL SUCH THAT: C C (-1)**(K-J) IF J .LE. K C P(X(J)) = C (-1)**(J+1+K) IF J .GT. K. C C A SUITABLE STARTING SEQUENCE FOR X CAN BE OBTAINED BY C CALLING SUBROUTINE SETUP. C N DEGREE OF EXTREMAL POLYNOMIAL. IF THIS IS TOO LARGE C OVERFLOW, UNDERFLOW, OR DIVISION BY ZERO ARE LIKELY C XINT TWO DIMENSIONAL ARRAY OF DIMENSION 2 BY NINT. ON INPUT, XINT C CONTAINS THE ENDPOINTS OF THE INTERVALS DEFINING S (SEE ABOVE C DEFINITION OF S). XINT MUST BE IN INCREASING ORDER: C XINT(1,1) .LE. XINT(2,1) .LT. XINT(1,2) .LE. XINT(2,2) .. C .LT. XINT(1,NINT) .LE. XINT(2,NINT). C XINT MUST BE SUCH THAT S CONTAINS AT LEAST N+1 POINTS. C NINT NUMBER OF INTERVALS IN S. C K INTEGER SUCH THAT X(K)=B AND X(K+1)=C. C IPRINT INTEGER WHICH ON INPUT HAS THE FOLLOWING MEANINGS C =-1 SUPPRESS ALL PRINTING. PZERO WILL NOT BE CALCULATED C =0 PZERO = P(0) WILL BE CALCULATED AND PRINTED C =1 PZERO AND THE ROOTS OF THE POLYNOMIAL WILL BE C CALCULATED AND PRINTED C =2 IN ADDITION, FOR EACH ITERATION THE CURRENT VALUE OF C THE ARRAY X WILL BE PRINTED C MAXIT MAXIMUM NUMBER OF ITERATIONS ALLOWED C EPS2 A TOLERANCE USED AS A STOPPING CRITERION FOR THE MAJOR C ITERATION. THE ROUTINE WILL STOP WHEN C MAX(ABS(XOLD(I)-X(I)),I=1,2,...,N+1) .LT. EPS2 C WHERE XOLD IS THE X ARRAY OF THE PREVIOUS ITERATION. C EPS3,EPS4 TOLERANCE PARAMETERS USED BY SUBROUTINE REGULA. THE C CALCULATED VALUE OF THE ROOT, XRT, WILL BE SUCH THAT XRT WILL C BE WITHIN EPS3 OF AN ACTUAL ROOT OR ABS(P(XRT)) WILL BE LESS C THAN EPS4. VALUES FOR EPS3 AND EPS4 NEED NOT BE PROVIDED IF C IPRINT IS EQUAL TO -1 OR 0. C XTILDE ARRAY OF DIMENSION N+1. ON OUTPUT, IF IPRINT IS EQUAL TO 1 OR C 2, XTILDE WILL CONTAIN THE RECIPROCALS OF THE ROOTS OF THE C EXTREMAL POLYNOMIAL. THE RECIPROCALS OF THE ROOTS IN THE C INTERVAL (X(1),X(N+1)) WILL BE STORED IN THE ELEMENTS C XTILDE(2),XTILDE(3),...,XTILDE(N). THE RECIPROCAL OF THE C REMAINING ROOT, IF IT EXISTS, WILL BE STORED IN XTILDE(1). C IF IT DOES NOT EXIST, XTILDE(1) WILL BE SET TO ZERO. C ALPHA ARRAY OF DIMENSION N+1 USED BY THIS ROUTINE TO STORE THE C VALUES (AS CALCULATED BY ALCALC) C C N+1 C ALPHA(I) = P(X(I)) / PRODUCT (X(I)-X(J)) C J=1 C J.NE.I C C WHERE P(X) IS THE VALUE OF THE POLYNOMIAL AT X. IN PARTICULAR, C C (-1)**(I-K) IF I .LE. K C P(X(I))= C (-1)**(I+1-K) IF I .GT. K C C ON RETURN, IF IPRINT IS UNEQUAL TO -1, ALPHA WILL CONTAIN C THESE VALUES FOR THE EXTREMAL POLYNOMIAL. C C PROD ARRAY OF DIMENSION N+1 USED BY THIS ROUTINE TO STORE C THE VALUES (AS CALCULATED BY ALCALC) C C N+1 C PROD(I) = PRODUCT (X(I) - X(J)) C J=1 C J.NE.I C C ON RETURN, IF IPRINT IS UNEQUAL TO -1, PROD WILL CONTAIN C THESE VALUES FOR THE EXTREMAL POLYNOMIAL. C C PZERO ON OUTPUT, IF IPRINT IS UNEQUAL TO -1, PZERO WILL CONTAIN C THE VALUE OF THE EXTREMAL POLYNOMIAL AT 0. C C IERR ON RETURN, IERR WILL BE THE SUM OF THE VALUES OF THE C FOLLOWING POSSIBLE USER INPUT ERRORS. IERR=0 IS THE C NORMAL RETURN. IF IERR IS NOT ZERO THE ROUTINE WILL C RETURN IMMEDIATELY WITHOUT CALCULATING THE EXTREMAL C POLYNOMIAL C C VALUE DESCRIPTION C ----- ----------- C 1 N .LT. 2 C 2 THE POINTS X ARE OUT OF ORDER (MUST BE INCREASING) C 4 NINT .LT. 2 C 8 THE POINTS OF XINT ARE OUT OF ORDER C 16 THE POINTS OF X ARE NOT IN S C 32 K .LT. 1, OR K .GT. N, OR X(K) .GE. 0, OR X(K+1) .LE. 0 C 64 EPS2 .LT. 0 C 128 EPS3 .LT. 0 AND IPRINT IS NOT -1 OR 0 C 256 EPS4 .LT. 0 AND IPRINT IS NOT -1 OR 0 C 512 THE POINT ZERO IS IN S C C OTHER SUBROUTINES C C SETUP NOT CALLED DIRECTLY OR INDIRECTLY BY EXTREM. MAY BE CALLED C BY THE USER PRIOR TO A CALL TO EXTREM TO OBTAIN A SUITABLE C STARTING SEQUENCE X. C ALCALC CALCULATES THE COEFFICIENTS ALPHA AND PROD OF A POLYNOMIAL C GIVEN X, M(=N+1), AND K. THE COEFFICIENTS ALPHA AND PROD ARE C USED IN SUBROUTINES EVAL, EVALI, AND DERIV TO EVALUATE THE C POLYNOMIAL AND ITS DERIVATIVE. C DERIV CALCULATES THE DERIVATIVE OF THE POLYNOMIAL AT THE POINT X(I) C FOR ANY I, 1 .LE. I .LE. N+1. C ERROR SUBROUTINE TO CHECK INPUT FOR POSSIBLE ERRORS. C EVAL EVALUATES THE POLYNOMIAL AT A GIVEN POINT. C EVALI EVALUATES (T**N)*P(1/T) FOR A GIVEN VALUE OF T. THIS IS USED C TO FIND THE ONE POSSIBLE ROOT OF P OUTSIDE OF THE INTERVAL C (X(1),X(N+1)). C FINDZR FINDS THE ROOTS OF THE POLYNOMIAL. C INTERP CALCULATES AN APPROXIMATION TO A MINIMA OF THE FUNCTION C -P(X(J))*P IN THE INTERVAL (X(J-1),X(J+1)) C REGULA SUBROUTINE CALLED BY FINDZR TO FIND A PARTICULAR ZERO USING C THE MODIFIED REGULA FALSI METHOD. C C C REMES ALGORITHM FOR THE EXTREMAL C C 1. CHOOSE (X(J),J=1,N+1) IN S, STRICTLY INCREASING, SUCH THAT X(K)=B, C X(K+1)=C FOR SOME K, WHERE B=MAX(T:T IN S, T.LE.0) AND C=MAX(T: C T IN S, T.GE.0) C C N+1 C 2. SET P(.) = SUM SIGN(L(J,0))*L(J,.) WITH C J=1 C C N+1 C L(J,T) = PRODUCT ((T-X(I))/(X(J)-X(I))) , ALL J C I=1 C I.NE.J C C 3. SET X(0) = MIN(T:T IN S), X(N+2) = MAX(T:T IN S) AND SET C C I X(J) FOR J=0,K,K+1,N+2 C I C U(J) = I THE FIRST OF THE POSSIBLY TWO MAXIMA OF C I P(X(J))*P(.) IN THE INTERSECTION OF S WITH THE C I CLOSED INTERVAL (X(J-1),X(J+1)), C I FOR J = 1,...,K-1,K+2,...,N+1 C C 4. CHOOSE THE VECTOR XTILDE FROM U AS FOLLOWS: C (A) IF P(X(0))*P(X(1)) .LT. -1, THEN XTILDE(I-1)=U(I), FOR C I=1,...,N+1, AND INCREASE K BY 1. C (B) IF P(X(N+2)*P(X(N+1) .LT. -1, THEN XTILDE(I)=U(I+1), FOR C I=1,...,N+1, AND DECREASE K BY 1. C (C) OTHERWISE, XTILDE(I) = U(I) , FOR I=1,...,N+1 C C 5. SET X(I)=XTILDE(I) , I=1,...,N+1 C C 6. ITERATE STEPS 2 THROUGH 5. C INTEGER K, N, NP1, ITER, IPRINT, IERR, MAXIT, ISET, NP2, KP1, * NINT, ISET0 REAL X(1), XTILDE(1), ALPHA(1), PROD(1), XINT(2,NINT) INTEGER I, J, JINT, JINT2 REAL EPS2, SIGN, SGSL, ERR, EPS3, PZERO, EPS4, DIFF, VAL, VAL2, * SLOPE CALL ERROR(X, N, XINT, NINT, K, IPRINT, EPS2, EPS3, EPS4, IERR) IF (IERR.NE.0) RETURN ISET0 = 0 NP1 = N + 1 NP2 = N + 2 ITER = 0 IF (N.LE.1) RETURN NOUT = 6 10 IF (IPRINT.EQ.2) WRITE (NOUT,99999) ITER, (X(I),I=1,NP1) 99999 FORMAT (7H0AFTER , I3, 28H ITERATIONS THE POINTS X ARE/(1X, * 10E13.5)) KP1 = K + 1 C GET THE INTERPOLATING POLYNOMIAL CALL ALCALC(X, NP1, K, ALPHA, PROD, XTILDE) SIGN = 1. IF (MOD(K,2).EQ.1) SIGN = -1. JINT = 1 ISET = 0 C IF ISET IS CHANGED THE RESULTING X'S MUST BE SHIFTED. IF ((X(1).EQ.XINT(1,1)) .OR. (ISET0.EQ.1)) GO TO 20 CALL EVAL(XINT(1,1), X, ALPHA, NP1, VAL) IF (VAL*SIGN.LE.1.) GO TO 20 ISET = -1 ISET0 = -1 20 DO 110 I=1,NP1 IF ((I.EQ.KP1) .OR. (I.EQ.K)) GO TO 100 C FIND JINT SO THAT X(I) IS IN THE CLOSED INTERVAL C (XINT(1,JINT),XINT(2,JINT)) 30 IF (X(I).LE.XINT(2,JINT)) GO TO 40 JINT = JINT + 1 GO TO 30 40 J = I - 1 CALL DERIV(X, I, PROD, ALPHA, NP1, SLOPE) SGSL = 1. IF (SLOPE*SIGN.GT.0.) SGSL = -1. IF (SGSL.EQ.1.) J = I + 1 IF (SLOPE.NE.0.) GO TO 50 XTILDE(I) = X(I) GO TO 110 50 SLOPE = -ABS(SLOPE) CALL INTERP(I, J, SIGN, X, ALPHA, NP1, SLOPE, SGSL, XINT, NINT, * XTILDE(I)) SIGN = -SIGN C C CHECK THAT XTILDE(I) IS IN S C JINT2 = JINT 60 IF (XTILDE(I).LE.XINT(2,JINT2)) GO TO 70 JINT2 = JINT2 + 1 GO TO 60 70 IF (XTILDE(I).GE.XINT(1,JINT2)) GO TO 80 JINT2 = JINT2 - 1 GO TO 70 80 IF (XTILDE(I).LE.XINT(2,JINT2)) GO TO 110 C HANDLE THE CASE WHEN XTILDE(I) IS NOT IN S CALL EVAL(XINT(1,JINT2+1), X, ALPHA, NP1, VAL) CALL EVAL(XINT(2,JINT2), X, ALPHA, NP1, VAL2) IF (VAL*SIGN.GT.VAL2*SIGN) GO TO 90 XTILDE(I) = XINT(2,JINT2) GO TO 110 90 XTILDE(I) = XINT(1,JINT2+1) GO TO 110 100 XTILDE(I) = X(I) SIGN = 1. 110 CONTINUE IF ((X(NP1).EQ.XINT(2,NINT)) .OR. (ISET0.EQ.(-1))) GO TO 120 CALL EVAL(XINT(2,NINT), X, ALPHA, NP1, VAL) IF (VAL*SIGN.GE.(-1.)) GO TO 120 ISET = 1 ISET0 = 1 C C SET UP X FOR NEXT ITERATION C 120 K = K - ISET ERR = 0. IF (ISET.NE.1) GO TO 130 ERR = XINT(2,NINT) - X(NP1) X(NP1) = XINT(2,NINT) 130 IF (ISET.NE.(-1)) GO TO 140 ERR = X(1) - XINT(1,1) X(1) = XINT(1,1) 140 DO 150 I=1,NP1 J = ISET + I IF ((J.EQ.0) .OR. (J.EQ.NP2)) GO TO 150 DIFF = ABS(X(I)-XTILDE(J)) IF (DIFF.GT.ERR) ERR = DIFF X(I) = XTILDE(J) 150 CONTINUE C CHECK IF DONE IF (ERR.LT.EPS2) GO TO 160 ITER = ITER + 1 IF (ITER.LT.MAXIT) GO TO 10 WRITE (NOUT,99998) MAXIT 99998 FORMAT (46H0MAXIMUM NUMBER OF ITERATIONS EXCEEDED, MAXIT=, I4) 160 IF (IPRINT.EQ.2) WRITE (NOUT,99997) (X(I),I=1,NP1) 99997 FORMAT (39H0AFTER FINAL ITERATION THE POINTS X ARE/(1X, 10E13.5)) IF (IPRINT.EQ.(-1)) RETURN C EVALUATE EXTREMAL POLYNOMIAL AT 0 . CALL ALCALC(X, NP1, K, ALPHA, PROD, XTILDE) CALL EVAL(0., X, ALPHA, NP1, PZERO) WRITE (NOUT,99996) PZERO 99996 FORMAT (7H0P(0) =, E20.8) IF (IPRINT.EQ.0) RETURN C FIND THE RECIPROCALS OF THE ZEROS OF THE EXTREMAL POLYNOMIAL CALL FINDZR(X, ALPHA, NP1, K, EPS3, EPS4, XTILDE) RETURN END SUBROUTINE ALCALC(X, M, K, ALPHA, PROD, WORK) ALC 10 C THIS SUBROUTINE CALCULATES THE COEFFICIENTS ALPHA AND PROD TO BE C USED TO EVALUATE THE LAGRANGE INTERPOLATING POLYNOMIAL AND ITS C DERIVATIVE. C C PARAMETERS C X ARRAY OF DIMENSION M CONTAINING THE SEQUENCE WHICH ALONG C WITH K AND M UNIQUELY DEFINES THE POLYNOMIAL. C M THE ORDER (DEGREE + 1) OF THE INTERPOLATING POLYNOMIAL. C K IS SUCH THAT ZERO IS IN THE INTERVAL (X(K),X(K+1)). C ALPHA AN ARRAY OF DIMENSION M. ON OUTPUT, ALPHA CONTAINS THE VALUES C M C ALPHA(I) = P(X(I)) / PRODUCT (X(I)-X(J)) C J=1 C J.NE.I C WHERE P(X) IS THE VALUE OF THE POLYNOMIAL AT X. IN PARTICULAR, C (-1)**(I-K) IF I .LE. K C P(X(I))= C (-1)**(I+1-K) IF I .GT. K C PROD AN ARRAY OF DIMENSION M. ON OUTPUT, PROD CONTAINS THE VALUES: C M C PROD(I) = PRODUCT (X(I) - X(J)) C J=1 C J.NE.I C C WORK AN ARRAY OF DIMENSION M USED AS A SCRATCH ARRAY. C INTEGER K, M, I, J, MM1 REAL X(M), ALPHA(M), PROD(M), WORK(M), XI, P, SIGN MM1 = M - 1 SIGN = 1. IF (MOD(K,2).EQ.0) SIGN = -1. DO 10 I=1,MM1 WORK(I) = X(I+1) 10 CONTINUE DO 30 J=1,M P = 1. XI = X(J) DO 20 I=1,MM1 P = P*(XI-WORK(I)) 20 CONTINUE PROD(J) = P WORK(J) = XI IF (J.EQ.K+1) SIGN = -SIGN ALPHA(J) = SIGN/P SIGN = -SIGN 30 CONTINUE RETURN END SUBROUTINE DERIV(X, I, PROD, ALPHA, M, D) DER 10 C THIS SUBROUTINE EVALUATES THE DERIVATIVE OF THE LAGRANGE INTER- C POLATING POLYNOMIAL AT THE I-TH POINT OF X. C C PARAMETERS C X ARRAY CONTAINING POINTS AT WHICH THE POLYNOMIAL IS KNOWN. C I THE DERIVATIVE IS TO BE EVALUATED AT X(I) C PROD AN ARRAY OF DIMENSION M WHICH CONTAINS COEFFICIENTS CALCU- C LATED BY SUBROUTINE ALCALC. C ALPHA AN ARRAY OF DIMENSION M WHICH CONTAINS COEFFICIENTS C CALCULATED BY SUBROUTINE ALCALC. C M ONE PLUS THE DEGREE OF THE POLYNOMIAL. C D ON OUTPUT, CONTAINS THE VALUE OF THE DERIVATIVE AT X(I). C C THE DERIVATIVE IS EVALUATED: C M C P'(X(I)) = PROD(I) * SUM ((ALPHA(J) + ALPHA(I))/(X(I) - X(J) C J=1 C J.NE.I C INTEGER I, M, J REAL X(M), PROD(M), D, XS, ALPHA(M), AS XS = X(I) AS = ALPHA(I) D = 0. DO 10 J=1,M IF (J.EQ.I) GO TO 10 D = D + (AS+ALPHA(J))/(XS-X(J)) 10 CONTINUE D = D*PROD(I) RETURN END SUBROUTINE ERROR(X, N, XINT, NINT, K, IPRINT, EPS2, EPS3, EPS4, ERR 10 * IERR) C THIS SUBROUTINE CHECKS THE PARAMETER FOR POSSIBLE USER INPUT C ERRORS. IERR RETURNED WILL BE THE SUM OF THE VALUE OF EACH C ERROR WHICH OCCURRED. C C VALUE DESCRIPTION C ----- ----------- C 1 N .LT. 2 C 2 THE POINTS X ARE OUT OF ORDER (MUST BE INCREASING) C 4 NINT .LT. 2 C 8 THE POINTS OF XINT ARE OUT OF ORDER C 16 THE POINTS OF X ARE NOT IN S C 32 K .LT. 1, OR K .GT. N, OR X(K) .GE. 0, OR X(K+1) .LE. 0 C 64 EPS2 .LT. 0 C 128 EPS3 .LT. 0 AND IPRINT IS NOT -1 OR 0 C 256 EPS4 .LT. 0 AND IPRINT IS NOT -1 OR 0 C 512 THE POINT ZERO IS IN S C INTEGER N, NINT, K, IPRINT, IERR, INT, I, NP1 REAL X(1), XINT(2,NINT), EPS2, EPS3, EPS4, XOLD IERR = 0 IF (N.GE.2) GO TO 10 C ERROR 1 IERR = IERR + 1 NOUT = 6 WRITE (NOUT,99999) 99999 FORMAT (38H ***ERROR*** N MUST BE GREATER THAN 1) IF (NINT.LT.2) GO TO 20 GO TO 140 10 IF (NINT.GE.2) GO TO 30 20 IERR = IERR + 4 WRITE (NOUT,99998) 99998 FORMAT (41H ***ERROR*** NINT MUST BE GREATER THAN 1) GO TO 140 30 IF ((K.LT.1) .OR. (K.GT.N)) GO TO 40 IF ((X(K).LT.0.) .AND. (X(K+1).GT.0.)) GO TO 50 40 IERR = IERR + 32 WRITE (NOUT,99997) 99997 FORMAT (43H ***ERROR*** K MUST BE BETWEEN 1 AND N, AND, 7H X(K) M, * 39HUST BE .LT. 0 AND X(K+1) MUST BE .GT. 0) GO TO 140 50 XOLD = X(1) - 1. INT = 1 IF (XINT(1,1).GT.XINT(2,1)) GO TO 110 NP1 = N + 1 DO 80 I=1,NP1 IF (X(I).LE.XOLD) GO TO 120 IF (X(I).LT.XINT(1,INT)) GO TO 130 IF (X(I).LE.XINT(2,INT)) GO TO 70 60 INT = INT + 1 IF (INT.GT.NINT) GO TO 130 IF (XINT(1,INT).LE.XINT(2,INT-1)) GO TO 110 IF (X(I).LT.XINT(1,INT)) GO TO 130 IF (X(I).GT.XINT(2,INT)) GO TO 60 70 XOLD = X(I) 80 CONTINUE C CHECK THAT ZERO IS NOT IN S DO 90 I=1,NINT IF (XINT(1,I).GT.0.) GO TO 140 IF (XINT(2,I).GE.0.) GO TO 100 90 CONTINUE 100 WRITE (NOUT,99996) 99996 FORMAT (35H ***ERROR*** ZERO MUST NOT BE IN S ) IERR = IERR + 512 GO TO 140 110 IERR = IERR + 8 WRITE (NOUT,99995) 99995 FORMAT (45H ***ERROR*** POINTS OF XINT ARE OUT OF ORDER ) GO TO 140 120 IERR = IERR + 2 WRITE (NOUT,99994) 99994 FORMAT (42H ***ERROR*** POINTS OF X ARE OUT OF ORDER ) GO TO 140 130 IERR = IERR + 16 WRITE (NOUT,99993) 99993 FORMAT (38H ***ERROR*** POINTS OF X ARE NOT IN S ) 140 IF (EPS2.GE.0.) GO TO 150 IERR = IERR + 64 WRITE (NOUT,99992) 99992 FORMAT (36H ***ERROR*** EPS2 IS LESS THAN ZERO ) 150 IF ((IPRINT.EQ.(-1)) .OR. (IPRINT.EQ.0)) RETURN IF (EPS3.GE.0.) GO TO 160 IERR = IERR + 128 WRITE (NOUT,99991) 99991 FORMAT (36H ***ERROR*** EPS3 IS LESS THAN ZERO ) 160 IF (EPS4.GE.0.) RETURN IERR = IERR + 256 WRITE (NOUT,99990) 99990 FORMAT (36H ***ERROR*** EPS4 IS LESS THAN ZERO ) RETURN END SUBROUTINE EVAL(T, X, ALPHA, M, P) EVA 10 C THIS SUBROUTINE EVALUATES THE VALUE OF THE LAGRANGE INTERPOLATING C POLYNOMIAL AT THE POINT T. C C PARAMETERS C T POINT AT WHICH POLYNOMIAL IS TO BE EVALUATED. C X ARRAY CONTAINING POINTS AT WHICH THE POLYNOMIAL IS KNOWN. C ALPHA ARRAY OF DIMENSION M WHICH HAS THE COEFFICIENTS CALCULATED C BY SUBROUTINE ALCALC. C M ONE PLUS THE DEGREE OF THE POLYNOMIAL. C P ON OUTPUT, THE VALUE OF THE POLYNOMIAL AT THE POINT T. C C THE POLYNOMIAL IS EVALUATED: C M M C P(T) = (PRODUCT(T-X(J))) * (SUM(ALPHA(J)/(T-X(J)))) C J=1 J=1 C INTEGER M, I REAL X(M), ALPHA(M), T, P, C, DIFF, SGN C = 1. P = 0. DO 10 I=1,M DIFF = (T-X(I)) C CHECK IF T EQUALS X(I) IN WHICH CASE P= 1 OR -1 IF (DIFF.EQ.0.) GO TO 20 C = C*DIFF P = P + ALPHA(I)/DIFF 10 CONTINUE P = P*C RETURN 20 SGN = 1. IF (MOD(M-I,2).EQ.1) SGN = -1. P = SIGN(1.,ALPHA(I)*SGN) RETURN END SUBROUTINE EVALI(T, X, ALPHA, M, P) EVA 10 C THIS SUBROUTINE EVALUATES (T**(M-1))*P(1/T) WHERE P(1/T) IS THE C LAGRANGE INTERPOLATING POLYNOMIAL EVALUATED AT 1/T C C PARAMETERS C T POINT AT WHICH THE ABOVE EXPRESSION IS TO BE EVALUATED. C X ARRAY CONTAINING POINTS AT WHICH THE POLYNOMIAL IS KNOWN. C ALPHA ARRAY OF DIMENSION M WHICH HAS THE COEFFICIENTS CALCULATED C BY SUBROUTINE ALCALC. C M ONE PLUS THE DEGREE OF THE POLYNOMIAL. C P ON OUTPUT, CONTAINS THE VALUE OF THE ABOVE EXPRESSION AT T. C C THE EXPRESSION IS EVALUATED: C M M C (T**(M-1))*P(1/T) =(PRODUCT(1-T*X(J)))*(SUM(ALPHA(J)/(1-T*X(J)))) C J=1 J=1 C INTEGER M, I REAL X(M), ALPHA(M), T, P, C, DIFF, SGN C = 1. P = 0. DO 10 I=1,M DIFF = (1.-T*X(I)) IF (DIFF.EQ.0.) GO TO 20 C IF DIFF EQUALS ZERO EVALUATE THE EXPRESSION DIFFERENTLY TO AVOID C DIVISION BY ZERO. C = C*DIFF P = P + ALPHA(I)/DIFF 10 CONTINUE P = P*C RETURN C P(1/T) IS EITHER +1 OR -1 20 SGN = 1. IF (MOD(M-I,2).EQ.1) SGN = -1. P = SIGN(1.,ALPHA(I)*SGN)*T**(M-1) RETURN END SUBROUTINE FINDZR(X, ALPHA, M, K, EPS3, EPS4, WORK) FIN 10 CALLS SUBROUTINES REGULA,EVAL,EVALI C C THIS SUBROUTINE FINDS THE ZEROS OF THE LAGRANGE INTERPOLATING POLY- C NOMIAL USING THE MODIFIED REGULA FALSI METHOD. THE RESULTS ARE C PRINTED OUT. C C PARAMETERS C X ARRAY OF DIMENSION M CONTAINING POINTS AT WHICH THE C VALUE OF THE POLYNOMIAL IS KNOWN. C ALPHA ARRAY OF DIMENSION M CONTAINING COEFFICIENTS CALCULATED BY C SUBROUTINE ALCALC. C M ORDER (=DEGREE +1) OF THE POLYNOMIAL. C K INTEGER SUCH THAT ZERO LIES IN THE INTERVAL X(K) TO X(K+1). C EPS3,EPS4 TOLERANCE PARAMETERS USED BY SUBROUTINE REGULA. THE C CALCULATED VALUE OF THE ROOT, XRT, WILL BE SUCH THAT C XRT WILL BE WITHIN EPS3 OF AN ACTUAL ROOT OR C ABS(P(XRT)) WILL BE LESS THAN EPS4. C WORK AN ARRAY OF DIMENSION M. ON OUTPUT, WORK CONTAINS C THE RECIPROCALS OF THE ROOTS. C C THERE WILL BE A ROOT BETWEEN EVERY TWO CONSECUTIVE POINTS OF X EXCEPT C BETWEEN X(K) AND X(K+1). THAT ACCOUNTS FOR M-2 ROOTS, THE RECIPROCALS C OF WHICH WILL BE PLACED IN WORK(2),WORK(3),...,WORK(M-1). THERE C MAY BE ONE MORE ROOT LYING OUTSIDE OF THE INTERVAL (X(1),X(M)). C IF THIS ROOT EXISTS THE RECIPROCAL WILL BE PLACE IN WORK(1),IF C NOT WORK(1) WILL BE SET TO ZERO. C INTEGER M, K, J, I, MM1 REAL X(M), WORK(M), ALPHA(M), SIGN, EPS3, EPS4 REAL A, B, FA, FB, SIGN1, SUM, XRT, ABSAL, ALMAX EXTERNAL EVAL, EVALI MM1 = M - 1 C C FIND THE M-2 ROOTS LYING IN THE INTERVAL (X(1),X(M)) C SIGN = 1. IF (MOD(K,2).EQ.1) SIGN = -1. SIGN1 = SIGN DO 10 J=1,MM1 IF (J.EQ.K) GO TO 10 A = X(J) B = X(J+1) FA = -SIGN FB = SIGN CALL REGULA(EVAL, A, B, FA, FB, X, ALPHA, M, EPS3, EPS4, XRT) SIGN = -SIGN I = J + 1 IF (J.GT.K) I = J WORK(I) = XRT 10 CONTINUE NOUT = 6 WRITE (NOUT,99999) (WORK(I),I=2,MM1) 99999 FORMAT (48H0THE ROOTS LYING IN THE INTERVAL (X(1),X(M)) ARE/(1X, * 5E20.7)) DO 20 J=2,M WORK(J) = 1./WORK(J) 20 CONTINUE C FIND ROOT OUTSIDE OF INTERVAL IF IT EXISTS SUM = 0. ALMAX = 0. DO 30 I=1,M ABSAL = ABS(ALPHA(I)) IF (ABSAL.GT.ALMAX) ALMAX = ABSAL SUM = SUM + ALPHA(I) 30 CONTINUE IF (ABS(SUM).GT.ALMAX*1.E-4) GO TO 50 40 WRITE (NOUT,99998) 99998 FORMAT (54H0THE ADDITIONAL ROOT CANNOT BE ACCURATELY CALCULATED O, * 13HR IS INFINITE) WORK(1) = 0. RETURN 50 IF (SUM*SIGN.GT.0.) GO TO 60 C THE ADDITIONAL ROOT IS TO THE LEFT OF X(1) A = 1./X(1) B = -1.E-4 FA = -SIGN1 CALL EVALI(B, X, ALPHA, M, FB) IF (FA*FB.GE.0.) GO TO 40 GO TO 70 C THE ADDITIONAL ROOT IS TO THE RIGHT OF X(M) 60 B = 1./X(M) A = 1.E-4 FB = -SIGN CALL EVALI(A, X, ALPHA, M, FA) IF (FA*FB.GE.0.) GO TO 40 70 CALL REGULA(EVALI, A, B, FA, FB, X, ALPHA, M, EPS3, EPS4, XRT) WRITE (NOUT,99997) XRT WORK(1) = XRT 99997 FORMAT (41H0THE RECIPROCAL OF THE ADDITIONAL ROOT IS, E20.7) RETURN END SUBROUTINE INTERP(I, J, SIGN, X, ALPHA, M, SLOPE, SGSL, XINT, INT 10 * NINT, XMIN) CALLS SUBROUTINE EVAL C C THIS SUBROUTINE FINDS AN APPROXIMATION TO XMIN WHERE XMIN IS SUCH C THAT: C C SIGN*P(XMIN)=MIN(SIGN*P(T):MIN(X(I),X(J)).LE.T.LE. MAX(X(I),X(J))) C C WHERE P(T) IS THE INTERPOLATING POLYNOMIAL AT THE POINT T. C C PARAMETERS C I,J INTEGERS SO THAT X(I) AND X(J) DEFINE THE ENDPOINTS OF THE C INTERVAL TO BE SEARCHED FOR XMIN. C SIGN REAL VALUE, EITHER +1. OR -1. SO THAT SIGN*P(X(I)) = -1. C X ARRAY OF DIMENSION M WHICH CONTAINS THE ABSCISSAE OF THE C POINTS USED FOR LAGRANGE INTERPOLATION. C ALPHA COEFFICIENTS FOR LAGRANGE POLYNOMIAL. C M ORDER (DEGREE PLUS ONE) OF INTERPOLATING POLYNOMIAL. C SLOPE -(ABS(P'(X(I)))) C SGSL HAS THE FOLLOWING VALUE C +1. IF X(I) .LT. X(J) C -1. IF X(I) .GT. X(J) C XINT MATRIX OF DIMENSION 2 BY NINT CONTAINING THE ENDPOINTS OF C THE INTERVALS DEFINING THE SET S. C NINT THE NUMBER OF INTERVALS IN S. C XMIN ON OUTPUT, CONTAINS THE CALCULATED VALUE OF XMIN. C C LET P(T) BE THE LAGRANGE INTERPOLATING POLYNOMIAL AT THE POINT T. C U, THE FIRST APPROXIMATION OF XMIN, IS THE POINT AT WHICH THE C PARABOLA DETERMINED BY P(X(I)), P(X(J)), AND P'(X(I)) TAKES ON ITS C MINIMUM. XMIN IS THEN TAKEN TO BE THE POINT AT WHICH THE PARABOLA C INTERPOLATING P(X(I)), P'(X(I)), AND P(U) TAKES ON ITS MINIMUM. C INTEGER I, J, M, NINT, IFIRST REAL SIGN, X(M), ALPHA(M), SLOPE, SGSL, XINT(2,NINT), XMIN REAL U, XMAX, F, XI, XT, DIFF2 XI = X(I) IFIRST = 1 IF (J.EQ.0) GO TO 10 IF (J.EQ.M+1) GO TO 20 U = ABS(XI-X(J)) XMAX = U F = 1. GO TO 40 C C HANDLE THE CASE WHEN X(I) IS AN ENDPOINT. 10 IF (XI.EQ.XINT(1,1)) GO TO 60 XT = XINT(1,1) GO TO 30 20 IF (XI.EQ.XINT(2,NINT)) GO TO 60 XT = XINT(2,NINT) C CALCULATE INTERPOLATING POLYNOMIAL AT XT 30 CALL EVAL(XT, X, ALPHA, M, F) F = F*SIGN U = ABS(XI-XT) XMAX = U C USE FOR FIRST STEP THE SLOPE AT X(I) AND THE POINT U. 40 DIFF2 = ((F+1.)/U-SLOPE)/U IF (DIFF2.LE.0.) GO TO 70 U = -SLOPE/DIFF2/2. IF (U.LE.XMAX) GO TO 50 XMIN = SGSL*XMAX + XI RETURN 50 XT = SGSL*U + XI CALL EVAL(XT, X, ALPHA, M, F) F = SIGN*F IF (IFIRST.NE.1) GO TO 70 IFIRST = 2 XMAX = U GO TO 40 60 XMIN = XI RETURN 70 XMIN = XI + SGSL*U RETURN END SUBROUTINE REGULA(FCT, A, B, FA, FB, X, ALPHA, M, EPS3, EPS4, XRT)REG 10 C THIS SUBROUTINE FINDS THE ZERO OF THE FUNCTION FCT WHICH IS C IN THE INTERVAL (A,B). THE MODIFIED REGULA FALSI METHOD IS USED. C C PARAMETERS C FCT THE NAME OF THE SUBROUTINE WHICH CALCULATES THE FUNCTION C FOR WHICH THE ROOT IS DESIRED. FCT IS EITHER EVAL (TO FIND C A ROOT OF THE POLYNOMIAL) OR EVALI (TO FIND THE RECIPROCAL C OF A ROOT OF THE POLYNOMIAL). FCT MUST BE DECLARED EXTERNAL C IN THE PROGRAM THAT CALLS REGULA. C A LEFT END POINT OF INTERVAL IN WHICH A ROOT LIES. C B RIGHT END POINT OF INTERVAL IN WHICH A ROOT LIES. C FA THE VALUE OF THE FUNCTION AT THE POINT A. C FB THE VALUE OF THE FUNCTION AT THE POINT B. C X ARRAY CONTAINING THE ABSCISSAE OF THE DATA POINTS. C ALPHA ARRAY CONTAINING COEFFICIENTS CALCULATED BY SUBROUTINE C ALCALC. C M ORDER OF THE POLYNOMIAL. C EPS3,EPS4 TOLERANCE PARAMETERS. THE CALCULATED VALUE OF XRT C WILL BE SUCH THAT XRT WILL BE WITHIN EPS3 OF AN C ACTUAL ROOT OR ABS(P(XRT)) WILL BE LESS THAN EPS4. C XRT ON OUTPUT, XRT CONTAINS THE CALCULATED VALUE OF THE ROOT. C INTEGER M REAL X(1), ALPHA(M), A, B, FA, FB, EPS3, EPS4, FW, XRT, FXRT EXTERNAL FCT FW = FA 10 XRT = (FB*A-FA*B)/(FB-FA) CALL FCT(XRT, X, ALPHA, M, FXRT) C CHECK IF THE FUNCTION IS CLOSE ENOUGH TO ZERO AT XRT. IF (ABS(FXRT).LE.EPS4) RETURN IF (SIGN(1.,FXRT)*FA.GT.0.) GO TO 20 B = XRT FB = FXRT IF (SIGN(1.,FXRT)*FW.GT.0.) FA = FA/2. GO TO 30 20 A = XRT FA = FXRT IF (SIGN(1.,FXRT)*FW.GT.0.) FB = FB/2. 30 FW = FXRT C CHECK IF A AND B ARE CLOSE ENOUGH TOGETHER TO STOP. IF (B-A.GT.EPS3) GO TO 10 XRT = (A+B)/2. RETURN END SUBROUTINE SETUP(N, XINT, NINT, NZERO, X, K, IERR, IWORK) SET 10 C THIS SUBROUTINE MAY BE CALLED BEFORE CALLING SUBROUTINE EXTREM TO C CALCULATE SUITABLE STARTING VALUES FOR X AND K. C C PARAMETERS C N THE DEGREE OF THE EXTREMAL POLYNOMIAL TO BE FOUND BY EXTREM. C XINT REAL ARRAY OF DIMENSION 2 BY NINT DEFINING THE SET S. S IS C THE UNION OF THE CLOSED INTERVALS (XINT(1,I),XINT(2,I)) FOR C I=1,2,...,NINT. XINT MUST BE INCREASING ORDER SUCH THAT: C XINT(1,I) .LE. XINT(2,I) I=1, 2, ..., NINT, AND C XINT(2,I) .LT. XINT(1,I+1) I=1, 2, ..., NINT-1. C NINT NUMBER OF INTERVALS IN S. C NZERO INTEGER SUCH THAT: C XINT(2,NZERO) .LE. 0 AND XINT(1,NZERO+1) .GE. 0 C X REAL ARRAY OF DIMENSION N+1. ON OUTPUT, X CONTAINS A STRICTLY C INCREASING SEQUENCE OF POINTS OF S. C K INTEGER VALUE, WHICH, ON OUTPUT, IS SUCH THAT: C X(K)=XINT(2,NZERO) AND X(K+1)=(1,NZERO+1) C IERR ON OUTPUT THIS ERROR FLAG CONTAINS ONE OF THE FOLLOWING: C 0 - NORMAL RETURN C 1 - INSUFFICIENT NUMBER OF POINTS IN S C 2 - ZERO NOT IN (XINT(2,NZERO),XINT(1,NZERO+1)) C 3 - XINT IS NOT IN INCREASING ORDER. C IWORK INTEGER WORK ARRAY OF DIMENSION AT LEAST NINT. C INTEGER N, NINT, NZERO, K, IWORK(NINT), I, NUMTOT, NM1, NP1, IW, J INTEGER NP1MNT, NUM, IERR, INDEX REAL XINT(2,NINT), X(1), TOTLEN, XLEN, DEL IF ((XINT(2,NZERO).GT.0.) .OR. (XINT(1,NZERO+1).LT.0.)) GO TO 130 IERR = 0 NP1 = N + 1 NM1 = N - 1 IF (NINT.GT.N) GO TO 90 C C CALCULATE THE TOTAL LENGTH OF THE DOMAIN. C TOTLEN = 0. DO 20 I=1,NINT IF (I.EQ.NINT) GO TO 10 IF (XINT(1,I+1).LE.XINT(2,I)) GO TO 140 10 XLEN = XINT(2,I) - XINT(1,I) IF (XLEN.LT.0.) GO TO 140 TOTLEN = TOTLEN + XLEN 20 CONTINUE IF (TOTLEN.EQ.0.) GO TO 120 C C CALCULATE APPROXIMATELY HOW MANY EXTREME POINTS SHOULD BE IN EACH C INTERVAL. C NUMTOT = 0 NP1MNT = N + 1 - NINT DO 30 I=1,NINT XLEN = XINT(2,I) - XINT(1,I) NUM = INT(FLOAT(NP1MNT)*(XLEN/TOTLEN)) NUMTOT = NUMTOT + NUM IWORK(I) = NUM + 1 IF (NUMTOT.LE.NM1) GO TO 30 IWORK(I) = IWORK(I) - (NUMTOT-NP1MNT) NUMTOT = NP1MNT 30 CONTINUE 40 IF (NUMTOT.EQ.NP1MNT) GO TO 60 DO 50 I=1,NINT IF (XINT(2,I).LE.XINT(1,I)) GO TO 50 IWORK(I) = IWORK(I) + 1 NUMTOT = NUMTOT + 1 IF (NUMTOT.EQ.NM1) GO TO 60 50 CONTINUE GO TO 40 60 NUMTOT = 0 C C IWORK(I) CONTAINS THE NUMBER OF POINTS OF X TO BE PLACED IN THE I-TH C INTERVAL OF S. DISTRIBUTE THESE POINTS EVENLY. C DO 80 I=1,NINT IW = IWORK(I) DEL = XINT(2,I) - XINT(1,I) IF (IW.GT.1) DEL = DEL/FLOAT(IW-1) DO 70 J=1,IW INDEX = NUMTOT + J X(INDEX) = XINT(1,I) + DEL*FLOAT(J-1) 70 CONTINUE NUMTOT = NUMTOT + IW IF (I.EQ.NZERO) K = NUMTOT 80 CONTINUE X(K) = XINT(2,NZERO) RETURN C C HANDLE THE CASE WHERE NINT IS GREATER THAN N. TRY TO DISTRIBUTE C THE POINTS AS UNIFORMLY AS POSSIBLE. C 90 NUMTOT = 0 DO 110 I=1,NINT IF (I.EQ.NZERO) GO TO 100 IF (I*NP1/NINT.LE.NUMTOT) GO TO 110 NUMTOT = NUMTOT + 1 X(NUMTOT) = XINT(1,I) GO TO 110 100 NUMTOT = NUMTOT + 2 IF (NUMTOT.GT.NP1) NUMTOT = NP1 K = NUMTOT - 1 X(K) = XINT(2,I) X(K+1) = XINT(1,I+1) 110 CONTINUE IF (NUMTOT.LE.N) X(N+1) = XINT(2,NINT) RETURN 120 NOUT = 6 WRITE (NOUT,99999) 99999 FORMAT (54H0ERROR IN SUBROUTINE SETUP, INSUFFICIENT NUMBER OF POI, * 18HNTS IN THE DOMAIN./35H EITHER INCREASE NINT OR DECREASE N) IERR = 1 RETURN 130 WRITE (NOUT,99998) 99998 FORMAT (54H0ERROR IN SUBROUTINE SETUP, ZERO IS NOT CONTAINED IN T, * 43HHE INTERVAL (XINT(2,NZERO),XINT(1,NZERO+1))) IERR = 2 RETURN 140 WRITE (NOUT,99997) 99997 FORMAT (54H0ERROR IN SUBROUTINE SETUP, XINT IS NOT IN INCREASING , * 5HORDER) IERR = 3 RETURN END C EXAMPLE PROGRAM FOR THE PACKAGE TO CONSTRUCT AN EXTREMAL POLYNOMIAL. MAN 10 C THIS EXAMPLE CONSISTS OF A SET S MADE UP OF 8 INTERVALS. TWO MAN 20 C EXTREMALS ARE TO BE CONSTRUCTED, THE FIRST OF DEGREE 4, THE MAN 30 C SECOND OF DEGREE 5. MAN 40 C SINCE S FOR THIS EXAMPLE IS SYMMETRIC, THE RESULTING EXTREMAL MAN 50 C POLYNOMIAL WILL BE AN EVEN FUNCTION. THEREFORE THE EXTREMAL MAN 60 C FOR N=4 WILL BE THE SAME AS FOR N=5. MAN 70 C MAN 80 C THE FOLLOWING IS A DESCRIPTION OF THE EXTREMAL POLYNOMIAL FOR THIS MAN 90 C CASE. MAN 100 C MAN 110 C THE POLYNOMIAL WILL HAVE ABSOLUTE VALUE OF ONE AT THESE SIX MAN 120 C POINTS OF S: MAN 130 C +/-.2, +/-SQRT(13/25), AND +/-1. MAN 140 C (SQRT(13/25) IS APPROXIMATELY .721110255) MAN 150 C MAN 160 C THE POLYNOMIAL WILL HAVE A VALUE OF 97/72 AT 0. MAN 170 C (97/72 IS APPROXIMATELY 1.347222222) MAN 180 C MAN 190 C THE POLYNOMIAL WILL HAVE ROOTS AT MAN 200 C +/-SQRT((13+6*SQRT(2))/25), AND +/-SQRT((13-6*SQRT(2))/25)) MAN 210 C (APPROXIMATELY +/-.927044365 AND +/-.424957345) MAN 220 C (1/SQRT((13+6*SQRT(2))/25) IS APPROXIMATELY 1.078697026) MAN 230 C MAN 240 C MAN 250 DIMENSION IWORK(8), XINT(2,8), X(6), XTILDE(6), ALPHA(6), PROD(6) MAN 260 C LEFT ENDPOINT , RIGHT ENDPOINT MAN 270 DATA XINT(1,1) /-1./, XINT(2,1) /-1./, XINT(1,2) /-.75/, MAN 280 * XINT(2,2) /-.35/, XINT(1,3) /-.3/, XINT(2,3) /-.3/, XINT(1,4) MAN 290 * /-.2/, XINT(2,4) /-.2/, XINT(1,5) /.2/, XINT(2,5) /.2/, MAN 300 * XINT(1,6) /.3/, XINT(2,6) /.3/, XINT(1,7) /.35/, XINT(2,7) MAN 310 * /.75/, XINT(1,8) /1./, XINT(2,8) /1./ MAN 320 DATA NZERO /4/, NINT /8/, N /4/, EPS2 /1.E-4/, EPS3 /1.E-5/, EPS4 MAN 330 * /1.E-5/ MAN 340 DATA MAXIT /30/, IPRINT /2/ MAN 350 C WRITE OUT THE DATA MAN 360 10 NOUT = 6 MAN 370 WRITE (NOUT,99999) N, EPS2, EPS3, EPS4 MAN 380 99999 FORMAT (10H1DEGREE = , I2, 9H, EPS2 = , E10.2, 9H, EPS3 = , MAN 390 * E10.2, 9H, EPS4 = , E10.2) MAN 400 WRITE (NOUT,99998) NINT, (XINT(1,I),XINT(2,I),I=1,NINT) MAN 410 99998 FORMAT (32H S IS THE UNION OF THE FOLLOWING, I3, 10H INTERVALS/ MAN 420 * 20H0 LEFT RIGHT /20H ENDPOINT ENDPOINT/(F6.2, 5X, F6.2))MAN 430 C CALL SUBROUTINE SETUP TO GET A STARTING SEQUENCE FOR X MAN 440 CALL SETUP(N, XINT, NINT, NZERO, X, K, IERR, IWORK) MAN 450 C CHECK FOR A NORMAL RETURN FROM SETUP MAN 460 IF (IERR.NE.0) GO TO 20 MAN 470 CALL EXTREM(X, N, XINT, NINT, K, IPRINT, MAXIT, EPS2, EPS3, EPS4, MAN 480 * XTILDE, ALPHA, PROD, PZERO, IERR) MAN 490 IF (IERR.NE.0) GO TO 20 MAN 500 IF (N.EQ.5) STOP MAN 510 N = 5 MAN 520 GO TO 10 MAN 530 20 WRITE (NOUT,99997) IERR MAN 540 99997 FORMAT (23H0ERROR IN DATA, IERR = , I2) MAN 550 STOP MAN 560 END MAN 570 C EXAMPLE PROGRAM TO DEMONSTRATE THE USE OF AN EXTREMAL POLYNOMIAL MAN 10 C TO OBTAIN ITERATION PARAMETERS FOR RICHARDSON ITERATION. MAN 20 C MAN 30 C THE PROBLEM TO BE SOLVED IS TO FIND X SUCH THAT A*X = RHS. MAN 40 C IN THIS EXAMPLE RHS IS A VECTOR OF ALL ONES, AND MAN 50 C A = (B**2 - SQRT(3)*I) WHERE B IS A 50 BY 50 TRIDIAGONAL MATRIX MAN 60 C WITH 2 ON THE DIAGONAL AND -1 ON THE OFF-DIAGONAL. MAN 70 C MAN 80 C THE EIGENVALUES OF A FOR THIS PROBLEM ARE MAN 90 C MAN 100 C (2. * (1. - COS(K*PI/51)))**2 - SQRT(3) K=1,...,50 MAN 110 C AND ARE THEREFORE ALL IN THE UNION OF THE TWO CLOSED INTERVALS MAN 120 C (-SQRT(3),-.2426) AND (.05094, 16-SQRT(3)). MAN 130 C MAN 140 C THIS EXAMPLE IS PROVIDED TO SHOW THE USE OF EXTREM, AND IS NOT MAN 150 C INTENDED AS A GENERAL SOLVER OF LINEAR SYSTEMS USING MAN 160 C RICHARDSON ITERATION. IN PARTICULAR, THE ORDERING OF MAN 170 C THE PARAMETERS HAS BEEN THE SUBJECT OF RESEARCH (E.G. SEE MAN 180 C V.I. LEBEDEV AND S.A. FINOGENOV (1971) ON THE ORDER OF CHOICE MAN 190 C OF THE ITERATION PARAMETERS IN THE CHEBYSHEV CYCLIC MAN 200 C ITERATION METHOD, ZHUR.VYCH.MAR. I FIZ. 11(2) 425-438. MAN 210 C AN ENGLISH TRANSLATION HAS APPEARED AS AN APPENDIX IN MAN 220 C R.S. ANDERSSEN AND G.H. GOLUB (1972) RHICHARDSONS NON- MAN 230 C STATIONARY MATRIX ITERATIVE PROCEDURE, REPORT MAN 240 C STAN-CS-72-304, COMP. SCI. DEPT., STANFORD U. ) MAN 250 C HOWEVER, IN THIS EXAMPLE THE ORDER OF THE PARAMETERS IS TAKEN MAN 260 C SO THAT IN EACH CYCLE THE I-TH ITERATION WILL USE PARAMETER MAN 270 C MOD(IJUMP*I,N)+1, WHERE N IS THE NUMBER OF PARAMETERS. (NOTE MAN 280 C THAT EACH PARAMETER WILL BE USED IF IJUMP AND N ARE RELATIVELY MAN 290 C PRIME. MAN 300 C MAN 310 C MAN 320 C SET UP THE PROBLEM FOR A MAXIMUM DEGREE OF 100, ALTHOUGH IF THE MAN 330 C DEGREE IS CHOSEN TOO LARGE OVERFLOW WILL OCCUR MAN 340 REAL XINT(2,2), X(100), ALPHA(100), PROD(100), XTILDE(100), R(50) MAN 350 INTEGER IWORK(2) MAN 360 INTEGER I, J, K, N, MAX, ICYCLE, INDEX, IPRINT, MAXIT, IJUMP, MAN 370 * NZERO MAN 380 INTEGER IERR, NINT MAN 390 REAL SUM, EPS2, EPS3, EPS4, SQR3, PZERO, RNORM, RMAX MAN 400 C ENTER THE SET IN WHICH THE SPECTRUM RESIDES MAN 410 DATA NINT /2/ MAN 420 DATA XINT(1,1) /-1.732051/, XINT(2,1) /-.2426/, XINT(1,2) MAN 430 * /.05094/, XINT(2,2) /14.267949/ MAN 440 C INITIALIZE PARAMETERS FOR EXTREM AND SETUP MAN 450 DATA IPRINT /1/, EPS2 /1.E-4/, EPS3 /1.E-5/, EPS4 /1.E-15/, MAXIT MAN 460 * /30/ MAN 470 DATA NZERO /1/ MAN 480 C GIVE VALUES FOR THE DEGREE, MAX NUMBER OF CYCLES FOR RICHARDSON MAN 490 C AND IJUMP. TO EXPERIMENT WITH THIS ROUTINE TAKE THE READ AND MAN 500 C WRITE STATEMENTS OUT OF COMMENTS BELOW TO READ IN N,MAX,IJUMP MAN 510 DATA N /10/, MAX /600/, IJUMP /3/ MAN 520 SQR3 = SQRT(3.) MAN 530 C WRITE(NOUT,1) MAN 540 C 1 FORMAT(27H INPUT N, MAX, IJUMP (3I5)) MAN 550 C READ(5,2) N, MAX, IJUMP MAN 560 C 2 FORMAT(3I5) MAN 570 CALL SETUP(N, XINT, NINT, NZERO, X, K, IERR, IWORK) MAN 580 WRITE (NOUT,99999) IERR MAN 590 99999 FORMAT (18H IERR FROM SETUP =, I3) MAN 600 CALL EXTREM(X, N, XINT, NINT, K, IPRINT, MAXIT, EPS2, EPS3, EPS4, MAN 610 * XTILDE, ALPHA, PROD, PZERO, IERR) MAN 620 WRITE (NOUT,99998) IERR MAN 630 99998 FORMAT (19H IERR FROM EXTREM =, I3) MAN 640 WRITE (NOUT,99997) (XTILDE(I),I=1,N) MAN 650 99997 FORMAT (30H0AFTER CALL TO EXTREM XTILDE= /(1X, 5E20.7)) MAN 660 C INITIALIZE X MAN 670 DO 10 I=1,50 MAN 680 X(I) = 0. MAN 690 10 CONTINUE MAN 700 DO 50 ICYCLE=1,MAX MAN 710 RMAX = 0. MAN 720 C CALCULATE RESIDUAL AND ITS L2 NORM MAN 730 DO 40 J=1,N MAN 740 SUM = 0. MAN 750 DO 20 I=3,48 MAN 760 R(I) = X(I-2) + X(I+2) - 4.*(X(I-1)+X(I+1)) + MAN 770 * (6.-SQR3)*X(I) - 1. MAN 780 SUM = SUM + R(I)**2 MAN 790 20 CONTINUE MAN 800 R(1) = (5.-SQR3)*X(1) - 4.*X(2) + X(3) - 1. MAN 810 R(2) = (6.-SQR3)*X(2) - 4.*(X(1)+X(3)) + X(4) - 1. MAN 820 R(49) = (6.-SQR3)*X(49) - 4.*(X(48)+X(50)) + X(47) - 1. MAN 830 R(50) = (5.-SQR3)*X(50) - 4.*X(49) + X(48) - 1. MAN 840 SUM = SUM + R(1)**2 + R(2)**2 + R(49)**2 + R(50)**2 MAN 850 RNORM = SQRT(SUM) MAN 860 IF (RNORM.GT.RMAX) RMAX = RNORM MAN 870 INDEX = MOD(IJUMP*J,N) + 1 MAN 880 DO 30 I=1,50 MAN 890 X(I) = X(I) - XTILDE(INDEX)*R(I) MAN 900 30 CONTINUE MAN 910 40 CONTINUE MAN 920 IF (MOD(ICYCLE,50).EQ.0) WRITE (NOUT,99996) ICYCLE, RMAX MAN 930 99996 FORMAT (13H CYCLE NUMBER, I4, 18H NORM OF RESIDUAL=, E20.7) MAN 940 50 CONTINUE MAN 950 STOP MAN 960 END MAN 970