C ALGORITHM 592, COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 1, C MAR., 1983, P. 98-116. SUBROUTINE RANGE(IAG, N, X, F, K, WK, LWK, AL, ALPHA, IL, RAN 10 * LEP, ETA, PSI, LOW, UP, OMEGA, IFAIL) C ******************************************************************** C PURPOSE * C ******************************************************************** C * C GIVEN VALUES OF A FUNCTION F(X) AT N DISTINCT POINTS * C X(1).LT.X(2),...,.LT.X(N) AND GIVEN A FINITE BOUND,AL, * C ON THE KTH. DERIVATIVE OF F(X), 1.LE.K.LE.N, * C THIS SUBROUTINE COMPUTES THE CLOSEST POSSIBLE BOUNDS * C ON F(ALPHA), WHERE ALPHA IS A SPECIFIED VALUE OF X. * C THE SUBROUTINE ALSO PROVIDES THE ESTIMATE, OMEGA, OF * C F(ALPHA) WHOSE ERROR HAS THE SMALLEST POSSIBLE BOUND. * C * C ******************************************************************** C C C **** I N P U T **** C C IAG IS AN INTEGER VARIABLE WHICH MUST BE SET TO THE VALUE +1 C AT THE FIRST CALL OF THE SUBROUTINE. THE SUBROUTINE MAY C BE RE-ENTERED WITH A DIFFERENT VALUE OF ALPHA. IN THIS C CASE IF THE VALUE OF IAG IS GREATER THAN 1, AND THE C REMAINING PARAMETERS ARE UNALTERED, THEN EXECUTION IS C MUCH FASTER.NOTE THAT THE CODE DOES NOT CHECK THAT THE C REMAINING PARAMETERS ARE UNALTERED. C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C N IS AN INTEGER VARIABLE WHICH MUST BE SET TO THE NUMBER C OF DATA POINTS X(I),I=1,...,N. RESTRICTION: N.GE.2 C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C X IS A REAL ARRAY OF LENGTH AT LEAST N WHICH MUST BE C SET TO THE VALUES OF THE DATA POINTS X(I),I=1,...,N. C RESTRICTION: THE DATA POINTS MUST BE DISTINCT AND C THEY MUST BE IN ASCENDING ORDER. C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C F IS A REAL ARRAY OF LENGTH AT LEAST N WHICH MUST BE C SET TO THE FUNCTION VALUES F(X(1)),...,F(X(N)). C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C K IS AN INTEGER VARIABLE WHICH MUST BE SET TO THE ORDER C OF THE DERIVATIVE OF F(X) FOR WHICH A FINITE BOUND IS C GIVEN. THE VALUE OF K SHOULD BE VERY MUCH SMALLER THAN C THE VALUE OF N. IN FACT WE RECOMMEND THAT ONLY IN C EXCEPTIONAL CIRCUMSTANCES AND THEN ONLY ON SOUND C NUMERICAL GROUNDS, SHOULD THE VALUE OF K BE GREATER C THAN 8. RESTRICTION: 1.LE.K.LE.N C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C WK IS A REAL ARRAY OF LENGTH AT LEAST: C C 5*N-2*K+1+(N-K)*MIN(K,N-K) C C WHICH IS USED AS WORKSPACE. C C LWK IS AN INTEGER VARIABLE WHICH MUST BE SET TO THE C LENGTH OF WK. C RESTRICTION: LWK.GE.5*N-2*K+1+(N-K)*MIN(K,N-K). C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C AL IS A REAL VARIABLE WHICH MUST BE SET TO THE VALUE C L, OF THE FINITE BOUND ON THE KTH. DERIVATIVE OF C F(X). C RESTRICTION: L MUST BE GREATER THAN THE LEAST VALUE C OF THE MAXIMUM ABSOLUTE VALUE OF THE KTH. DERIVATIVE C OF F(X) THAT IS CONSISTENT WITH THE GIVEN FUNCTION C VALUES F(X(1)),...,F(X(N)). IN PARTICULAR L MUST C SATISFY THE INEQUALITIES C C L .GT. 0 C AND C L .GE. FACTORIAL(K)*ABS(F(X(I),...,X(I+K)) C I=1,...,N-K, C C WHERE F(X(I),...,X(I+K)) DENOTES THE KTH. DIVIDED C DIFFERENCE OF F(X) BASED ON THE POINTS X(I),...,X(I+K). C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C ALPHA IS A REAL VARIABLE WHICH MUST BE SET TO THE VALUE C OF THE ARGUMENT X AT WHICH THE RANGE OF POSSIBLE C VALUES OF F(X) IS COMPUTED. C RESTRICTION: X(1).LE.ALPHA C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C IL IS AN INTEGER ARRAY OF LENGTH AT LEAST N. IT IS USED C AS WORKSPACE. C C LEP IS AN INTEGER VARIABLE WHICH MUST BE SET TO THE LESSER C LENGTH OF ARRAYS ETA AND PSI. RESTRICTION: LEP.GE.N-K. C THIS ARGUMENT IS NOT ALTERED BY THE SUBROUTINE. C C C **** O U T P U T **** C C ETA IS A REAL ARRAY OF LENGTH AT LEAST N-K. ON EXIT C FROM THE SUBROUTINE ETA CONTAINS THE KNOTS C OF THE PERFECT SPLINE U(X). C C PSI IS A REAL ARRAY OF LENGTH AT LEAST N-K. ON EXIT C FROM THE SUBROUTINE PSI CONTAINS THE KNOTS OF THE C PERFECT SPLINE L(X). C C LOW IS A REAL VARIABLE. ON EXIT FROM THE SUBROUTINE C LOW IS SET TO THE GREATEST LOWER BOUND OF F(ALPHA). C C UP IS A REAL VARIABLE. ON EXIT FROM THE SUBROUTINE C UP IS SET TO THE LEAST UPPER BOUND OF F(ALPHA). C C OMEGA IS A REAL VARIABLE. ON EXIT FROM THE SUBROUTINE C OMEGA IS SET TO THE OPTIMAL ESTIMATE OF F(ALPHA). C THE SMALLEST VALUE OF THE ERROR OF THIS ESTIMATE C OF F(ALPHA) IS ZERO, AND THE MAXIMUM VALUE WHICH C IT CAN ATTAIN IS THE QUANTITY: 0.5*ABS(UP-LOW). C C IFAIL IS AN ERROR RETURN FLAG. ON EXIT FROM THE SUBROUTINE C IT HAS ONE OF THE FOLLOWING VALUES: C C 0 SUCCESSFUL ENTRY C 1 N .LT. 2 C 2 K .LT. 1 OR K .GT. N C 3 X(I) .GE. X(I+1) FOR SOME I C 4 L .LT. FACTORIAL(K)*ABS(F(X(I),...,X(I+K)) C FOR SOME I C 5 MORE THAN IMAX ITERATIONS NEEDED TO CALCULATE C THE KNOTS ETA AND/OR PSI. C 6 THE METHOD USED TO CALCULATE THE KNOTS ETA, C AND PSI, FAILED. THIS IS USUALLY BECAUSE L C IS TOO SMALL OR THE VALUE OF K IS TOO LARGE. C 7 L .LE. 0 C 8 THE ARRAY WK IS TOO SMALL C 9 THE ARRAYS ETA AND PSI ARE TOO SMALL C 10 ALPHA .LT. X(1) C C C **** ADDITIONAL ROUTINES **** C C THE FOLLOWING ADDITIONAL ROUTINES ARE SUPPLIED WITH RANGE: C C SUBROUTINE KNOTS C SUBROUTINE RESID C SUBROUTINE JAC C C SINCE THESE ROUTINES ARE CALLED FROM WITHIN RANGE THE USER C SHOULD ENSURE THAT THERE ARE NO POTENTIAL PROBLEMS DUE TO C NAME CONFLICTS. C C C **** QUALITY ASSURANCE AND SOFTWARE STANDARD **** C C THE SUBROUTINES THAT COMPRISE THIS PACKAGE C HAVE BEEN WRITTEN TO CONFORM TO THE FORTRAN IV C ANSI STANDARD 1966, AND THEY HAVE BEEN VERIFIED C USING THE BELL TELEPHONE LABORATORIES FORTRAN C VERIFIER: PFORT. C THE SUBROUTINES HAVE BEEN EXTENSIVELY TESTED ON C A VARIETY OF TEST PROBLEMS, AND THEY HAVE BEEN C ANALYSED FOR ERRORS USING THE DAVE SYSTEM FROM C THE UNIVERSITY OF COLORADO. C TO MAKE THE CODE EASY TO READ THE SUBROUTINES C HAVE BEEN REFORMATTED USING POLISH. C C C C **** P.W.GAFFNEY DECEMBER 30. 1981 **** C C ******************************************************************** C REAL X(N), F(N), WK(LWK), ETA(LEP), PSI(LEP), LOW REAL LAMBDA INTEGER IL(N), HALFK C C ************************************************************** C N.B. IN THE FOLLOWING DATA STATEMENT * C IMAX IS A LIMIT IMPOSED ON THE NUMBER OF NEWTON * C ITERATIONS REQUIRED TO COMPUTE THE KNOTS ETA AND * C PSI. * C DEFAULT: IMAX=1000 * C NOUT IS THE LOGICAL UNIT NUMBER FOR THE OUTPUT UNIT * C WHERE ERROR MESSAGES ARE PRINTED. * C DEFAULT: NOUT = 6 * C TO SUPPRESS PRINTING SET NOUT .LE. 0 * C * DATA IMAX, NOUT /1000,6/ C ************************************************************** C DATA ZERO, POINT5, ONE /0.0,0.5,1.0/ IFAIL = 0 NMK = N - K NK = MIN0(K,NMK) LAJ = NMK*NK KM1 = K - 1 KP1 = K + 1 KP2 = K + 2 NP1 = N + 1 NP2 = N + 2 C*************************************************************** C CHECK THAT THE INPUT PARAMETERS ARE SENSIBLE * C*************************************************************** IF (ALPHA.LT.X(1)) GO TO 400 IF (N.LT.2) GO TO 310 IF (LWK.LT.5*N-2*K+1+LAJ) GO TO 380 IF (LEP.LT.NMK) GO TO 390 IF (K.LT.1 .OR. K.GT.N) GO TO 320 IF (AL.LE.ZERO) GO TO 370 NM1 = N - 1 DO 10 I=1,NM1 IF (X(I).LT.X(I+1)) GO TO 10 GO TO 330 10 CONTINUE C C COMPUTE FACTORIAL K-1. NOTE THAT IF K HAS A VALUE MUCH LARGER C THAN 8 THEN THERE IS THE POTENTIAL FOR OVERFLOW IN THE FOLLOWING C CALCULATION. C FACT = ONE IF (KM1.LT.2) GO TO 30 DO 20 I=2,KM1 FACT = FACT*FLOAT(I) 20 CONTINUE C END OF COMPUTING FACTORIAL K-1. C 30 IF (IAG.GT.1) GO TO 80 IF (K.EQ.N) GO TO 80 C C COMPUTE THE KTH. DIVIDED DIFFERENCES F(X(I),...,X(I+K)) C AND STORE THEM IN WK(K+1+I),I=1,...,N-K C DO 40 I=1,N KI = KP1 + I WK(KI) = F(I) 40 CONTINUE DO 60 J=1,K C J DENOTES THE JTH DIVIDED DIFFERENCE OF F. NMJ = N - J DO 50 I=1,NMJ C WK(KP1+I) IS THE DIVIDED DIFFERENCE F(X(I),...,X(I+J-1)). KI = KP1 + I S1 = WK(KI+1) - WK(KI) IPJ = I + J S1 = S1/(X(IPJ)-X(I)) WK(KI) = S1 50 CONTINUE 60 CONTINUE C C END OF COMPUTING THE DIVIDED DIFFERENCES. C C COMPUTE THE COMPONENTS OF CSTA, NAMELY: C FACTORIAL(K-1)*F(X(I),...,X(I+K))/AL, I=1,...,N-K, C IN WK(K+1+I),I=1,...,NMK, AND NOTE THE LARGEST COMPONENT C IN THE VARIABLE WKM. C WKM = ZERO DO 70 I=1,NMK KI = KP1 + I WK(KI) = FACT*WK(KI)/AL WKM = AMAX1(ABS(WK(KI)),WKM) 70 CONTINUE C C END OF COMPUTING THE COMPONENTS OF CSTA. C C CHECK THAT THE VALUE OF AL IS LARGE ENOUGH C IF (WKM.GT.ONE/FLOAT(K)) GO TO 340 C*************************************************************** C END OF CHECKING THE INPUT PARAMETERS. * C*************************************************************** C * C*************************************************************** C COMPUTE MACHINE EPSILON * C*************************************************************** 80 EPS = ONE 90 EPS = POINT5*EPS EPSP1 = EPS + ONE IF (EPSP1.GT.ONE) GO TO 90 EPS = EPS + EPS C*************************************************************** C END OF COMPUTING MACHINE EPSILON. * C*************************************************************** C C*************************************************************** C COMPUTE JINT SUCH THAT X(JINT).LE.ALPHA.LT.X(JINT+1) * C NOTE THAT THE FOLLOWING CODE MAKES USE OF THE FACT THAT * C IN FORTRAN THE RESULT OF INTEGER DIVISION IS ROUNDED * C DOWN TO THE NEAREST INTEGER. * C*************************************************************** ILEFT = 1 IRITE = NP1 100 MIDDLE = (ILEFT+IRITE)/2 IF (ALPHA.LT.X(MIDDLE)) IRITE = MIDDLE IF (ALPHA.GE.X(MIDDLE)) ILEFT = MIDDLE IF (IRITE.GT.ILEFT+1) GO TO 100 JINT = ILEFT C*************************************************************** C END OF COMPUTING JINT. * C*************************************************************** C IF (K.GT.1) GO TO 140 C C*************************************************************** C SET THE SOLUTION FOR K=1 AND RETURN * C*************************************************************** C C IF RE-ENTERED,(WITH ALPHA AND IAG THE ONLY CHANGED C PARAMETERS), THEN DO NOT RE-COMPUTE THE KNOTS C IF (IAG.GT.1) GO TO 120 DO 110 I=1,NMK IP1 = I + 1 XP = X(I) + POINT5*(X(IP1)-X(I)) FP = F(IP1) - F(I) ETA(I) = XP + POINT5*(FLOAT((-1)**IP1)*FP)/AL PSI(I) = XP + POINT5*(FLOAT((-1)**I)*FP)/AL 110 CONTINUE 120 I = JINT J = JINT IF (ALPHA.GE.X(N)) GO TO 130 IF (ALPHA.GE.ETA(JINT)) I = JINT + 1 IF (ALPHA.GE.PSI(JINT)) J = JINT + 1 130 TEMP1 = FLOAT((-1)**(I-1))*AL*(ALPHA-X(I)) + F(I) TEMP2 = FLOAT((-1)**(J))*AL*(ALPHA-X(J)) + F(J) LOW = AMIN1(TEMP1,TEMP2) UP = AMAX1(TEMP1,TEMP2) OMEGA = POINT5*(LOW+UP) C*************************************************************** C RETURN TO CALLING PROGRAM * C*************************************************************** C GO TO 410 C C*************************************************************** C COMPUTE THE CLOSEST POSSIBLE BOUNDS ON F(ALPHA) FOR K.KE.2 * C*************************************************************** C C IF RE-ENTERED,(WITH ALPHA AND IAG THE ONLY CHANGED C PARAMETERS), THEN DO NOT RE-COMPUTE THE KNOTS C 140 IF (IAG.GT.1) GO TO 160 C IF (K.EQ.N) GO TO 160 C C*************************************************************** C COMPUTE THE KNOTS OF THE PERFECT SPLINES * C*************************************************************** IWK = 4*NMK + (2*K) INCP1 = LWK - LAJ + 1 C C FIRST COMPUTE THE KNOTS ETA(1),...,ETA(N-K) CALL KNOTS(N, X, K, WK(KP2), WK(INCP1), NK, WK(NP2), IWK, * IL, NMK, ETA, IMAX, IFAIL, EPS) C C THEN(IF IFAIL=0) COMPUTE THE KNOTS PSI(1),...,PSI(N-K) IF (IFAIL.EQ.5) GO TO 350 IF (IFAIL.EQ.6) GO TO 360 DO 150 I=1,NMK KI = KP1 + I WK(KI) = -WK(KI) 150 CONTINUE CALL KNOTS(N, X, K, WK(KP2), WK(INCP1), NK, WK(NP2), IWK, * IL, NMK, PSI, IMAX, IFAIL, EPS) IF (IFAIL.EQ.5) GO TO 350 IF (IFAIL.EQ.6) GO TO 360 C*************************************************************** C END OF COMPUTING THE KNOTS. * C*************************************************************** C C*************************************************************** C SET THE KNOTS OF THE AUGMENTED B-SPLINE MALPHA(X),BY * C CHOOSING K DATA POINTS CLOSEST TO ALPHA AND ARRANGING * C THE RESULTING K+1 POINTS IN ASCENDING ORDER. * C STORE THE KNOTS IN WK(1),...,WK(K+1), AND STORE * C THE SUBSCRIPT OF THE DATA POINT WHICH FORMS THE FIRST * C KNOT OF MALPHA(X) IN THE VARIABLE IBAR. * C*************************************************************** 160 HALFK = K/2 IF (JINT.LE.HALFK) GO TO 190 IF (JINT.GE.NMK+HALFK) GO TO 220 J1 = JINT - HALFK DO 170 I=1,HALFK II = J1 + I WK(I) = X(II) 170 CONTINUE WK(HALFK+1) = ALPHA J2 = HALFK + 2 DO 180 I=J2,KP1 II = J1 - 1 + I WK(I) = X(II) 180 CONTINUE IBAR = J1 + 1 GO TO 250 190 DO 200 I=1,JINT WK(I) = X(I) 200 CONTINUE WK(JINT+1) = ALPHA J2 = JINT + 2 DO 210 I=J2,KP1 WK(I) = X(I-1) 210 CONTINUE IBAR = 1 GO TO 250 220 IEND = JINT - NMK DO 230 I=1,IEND II = NMK + I WK(I) = X(II) 230 CONTINUE WK(IEND+1) = ALPHA J2 = IEND + 2 IBAR = NMK + 1 IF (J2.GT.KP1) GO TO 250 DO 240 I=J2,KP1 II = NMK - 1 + I WK(I) = X(II) 240 CONTINUE C*************************************************************** C END OF SETTING THE KNOTS OF MALPHA(X). * C*************************************************************** C C*************************************************************** C COMPUTE THE QUANTITIES CUP AND CLOW * C*************************************************************** C C SET CUP AND CLOW WHEN K=N 250 IF (K.LT.N) GO TO 260 CUP = AL/FLOAT(K) CLOW = -CUP GO TO 270 C C CALCULATE CUP AND CLOW WHEN K.LT.N C FIRST COMPUTE THE INTEGRAL OF MALPHA*(THE KTH.DERIV. OF U(X)/AL) 260 I2 = NP2 + K CALL RESID(-1, 1, 1, KP1, WK(1), NMK, ETA, WK(KP2), K, * WK(NP2), WK(I2), NP2, WK, NP1, IL(KP1), SS, ICONT, * LAMBDA) TEMP = AL*WK(NP2) C C THEN COMPUTE THE INTEGRAL OF MALPHA*(THE KTH.DERIV. OF L(X)/AL) CALL RESID(-1, 1, 1, KP1, WK(1), NMK, PSI, WK(KP2), K, * WK(NP2), WK(I2), NP2, WK, NP1, IL(KP1), SS, ICONT, * LAMBDA) CUP = AMAX1(TEMP,-AL*WK(NP2)) CLOW = AMIN1(TEMP,-AL*WK(NP2)) C*************************************************************** C END OF COMPUTING CUP AND CLOW. * C*************************************************************** C C*************************************************************** C USING A MODIFIED NEVILLE ALGORITHM COMPUTE THE VALUE, * C AT X=ALPHA, OF THE POLYNOMIAL OF DEGREE K-1 WHICH * C PASSES THROUGH THE K FUNCTION VALUES F(X(IBAR)),..., * C F(X(IBAR+K-1)), AND STORE THE RESULT IN THE VARIABLE * C PALPHA. * C*************************************************************** C 270 WK(1) = F(IBAR) DO 290 I=2,K IBI = IBAR + I - 1 WK(I) = F(IBI) IM1 = I - 1 DO 280 IJ=1,IM1 J = I - IJ IBJ = IBAR + J - 1 WK(J) = WK(J+1) + (WK(J+1)-WK(J))*(ALPHA-X(IBI)) * /(X(IBI)-X(IBJ)) 280 CONTINUE 290 CONTINUE PALPHA = WK(1) C*************************************************************** C END OF EVALUATING THE INTERPOLATING POLYNOMIAL. * C*************************************************************** C C*************************************************************** C COMPUTE THE RANGE OF POSSIBLE VALUES OF F(ALPHA) * C AND COMPUTE THE OPTIMAL ESTIMATE OMEGA * C*************************************************************** C PROD = ONE DO 300 I=1,K LI = IBAR + I - 1 PROD = PROD*(ALPHA-X(LI)) 300 CONTINUE TEMP1 = PALPHA + PROD*CLOW/FACT TEMP2 = PALPHA + PROD*CUP/FACT LOW = AMIN1(TEMP1,TEMP2) UP = AMAX1(TEMP1,TEMP2) OMEGA = POINT5*(LOW+UP) C*************************************************************** C END OF COMPUTING THE CLOSEST POSSIBLE BOUNDS ON F(ALPHA)* C*************************************************************** GO TO 410 C*************************************************************** C DIAGNOSTIC PRINTING * C*************************************************************** 310 IFAIL = 1 IF (NOUT.GT.0) WRITE (NOUT,99999) GO TO 410 320 IFAIL = 2 IF (NOUT.GT.0) WRITE (NOUT,99998) GO TO 410 330 IFAIL = 3 IF (NOUT.GT.0) WRITE (NOUT,99997) GO TO 410 340 IFAIL = 4 WKM = WKM*FLOAT(K)*AL IF (NOUT.GT.0) WRITE (NOUT,99996) WKM GO TO 410 350 IF (NOUT.GT.0) WRITE (NOUT,99995) GO TO 410 360 IF (NOUT.GT.0) WRITE (NOUT,99993) GO TO 410 370 IFAIL = 7 IF (NOUT.GT.0) WRITE (NOUT,99994) GO TO 410 380 IFAIL = 8 IF (NOUT.GT.0) WRITE (NOUT,99992) GO TO 410 390 IFAIL = 9 IF (NOUT.GT.0) WRITE (NOUT,99991) GO TO 410 400 IFAIL = 10 IF (NOUT.GT.0) WRITE (NOUT,99990) C END OF PRINTING. C*************************************************************** 410 RETURN C*************************************************************** 99999 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HN.LT.2 **, * 3H***) 99998 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HK IS NOT , * 6HIN THE, 28H INTERVAL 1.LE.K.LE.N *****) 99997 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HTHE DATA , * 6HPOINTS, 34H ARE NOT IN ASCENDING ORDER. *****) 99996 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HTHE VALUE, * 6H OF AL, 17H MUST BE AT LEAST, 1PE14.6, 6H *****) 99995 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HMETHOD FA, * 6HILED T, 36HO CONVERGE IN IMAX ITERATIONS. *****) 99994 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HTHE VALUE, * 6H OF AL, 19H IS .LE. ZERO *****) 99993 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HTHE CONTI, * 6HNUATIO, 20HN METHOD HAS FAILED.//1X, 10HTHIS IS US, * 10HUALLY BECA, 22HUSE AL IS TOO SMALL OR//1X, * 18HTHE VALUE OF K IS , 16HTOO LARGE. *****) 99992 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9H THE ARRA, * 6HY WK I, 17HS TOO SMALL *****) 99991 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9H THE ARRA, * 6HYS ETA, 28H AND PSI ARE TOO SMALL *****) 99990 FORMAT (1X, 24H***** MESSAGE FROM RANGE, 3X, 9HALPHA.LT., * 10HX(1) *****) END SUBROUTINE KNOTS(N, X, K, CSTA, AJ, NK, WK, ISW, L, NMK, KNO 10 * ETA, IMAX, IFAIL, EPS) C *************************************************************** C * PURPOSE * C *************************************************************** C * * C * GIVEN THE DATA POINTS: X(1).LT.X(2).,...,.LT.X(N) * C * AND THE VECTOR: (CSTA(1),...,CSTA(N-K)), WHERE 1.LE.K.LE.N-1* C * THIS SUBROUTINE COMPUTES THE KNOTS OF THE PERFECT SPLINE U * C * * C *************************************************************** C REAL X(N), CSTA(NMK), AJ(NMK,NK), WK(ISW), ETA(NMK), * LASTSS, LAMBDA INTEGER L(NMK) DATA ZERO, POINT5, ONE, TWO, FIVE, TON * /0.0,0.5,1.0,2.0,5.0,100.0/ IFAIL = 0 LAMBDA = ONE N2 = 2*NMK N3 = 3*NMK N1P1 = NMK + 1 N3P1 = N3 + 1 I1 = 4*NMK + 1 I2 = I1 + K EPS1 = TON*FLOAT(NMK*NMK)*EPS*EPS C*************************************************************** C SET THE INITIAL APPROXIMATION * C*************************************************************** KP1 = K + 1 DO 20 I=1,NMK SUM = ZERO DO 10 J=1,KP1 IPJM1 = I + J - 1 SUM = SUM + X(IPJM1) 10 CONTINUE ETA(I) = SUM/FLOAT(KP1) 20 CONTINUE C END OF INITIALISATION. C*************************************************************** C C THE VALUES OF THE INTEGER VARIABLES: ICONT,IUP,IARG AND IA C C ICONT HAS THE VALUE -1 WHEN A CONTINUATION STEP IS IN C PROGRESS. OTHERWISE IT HAS THE VALUE +1. C C IUP IS SET TO +1 WHEN ETA IS UPDATED AND -1 OTHERWISE. C C IARG IS USED IN RESID. THE RESIDUALS ARE EVALUATED AT ETA C IF IARG=+1 AND AT ETA+(NEWTON CORRECTION) IF IARG=-1 C C IA IS DEFINED IN RESID. IN THIS SUBROUTINE C IT IS SET TO THE VALUE +1. C C ICONT = 1 IUP = -1 IA = 1 C C*************************************************************** C OUTER LOOP FOR CONTINUATION METHOD * C*************************************************************** DO 230 JCONT=1,IMAX IARG = 1 C C*************************************************************** C BEGIN NEWTON ITERATION * C*************************************************************** DO 160 ITER=JCONT,IMAX C C COMPUTE THE RESIDUALS IN WK(NMK+1),...,WK(2*NMK), AND C COMPUTE THE SUM OF SQUARES SS CALL RESID(IA, IARG, NMK, N, X, NMK, ETA, CSTA, * K, WK(I1), WK(I2), N3, WK, NMK, L, SS, * ICONT, LAMBDA) IF (ITER.EQ.JCONT) GO TO 50 C C IF SS HAS DECREASED THEN UPDATE THE APPROXIMATION IF (SS.LT.LASTSS) GO TO 30 C C IF SS HAS INCREASED DUE TO ROUNDING ERROR FINISH ITERATION IF (SS.LE.LASTSS+FIVE*EPS*EPS) GO TO 140 C C ELSE SET LAMBDA AND BEGIN CONTINUATION STEP LAMBDA = POINT5*LAMBDA IF (LAMBDA.LE.EPS) GO TO 250 GO TO 170 C C UPDATE THE APPROXIMATION 30 IUP = 1 DO 40 I=1,NMK ETA(I) = ETA(I) + WK(I) 40 CONTINUE C END OF UPDATE. C C IF SS IS LESS THAN OR EQUAL TO EPS1 FINISH ITERATION IF (SS.LE.EPS1) GO TO 140 C C STORE THE SUM OF SQUARES 50 LASTSS = SS C C COMPUTE THE NEWTON CORRECTION IN WK(1),...,WK(NMK) CALL JAC(N, X, NMK, ETA, L, K, WK(I1), NK, AJ, * WK(1), WK(N1P1)) IARG = -1 C C C COMPUTE THE SQUARE OF THE L2-NORM OF THE CORRECTION AND OF ETA C AND THEN TEST FOR CONVERGENCE S1 = ZERO S2 = ZERO DO 60 I=1,NMK S1 = S1 + WK(I)*WK(I) S2 = S2 + ETA(I)*ETA(I) 60 CONTINUE IF (S1.LE.EPS*EPS*S2) GO TO 140 C C TEST THAT THE COMPONENTS OF ETA+(NEWTON CORRECTION) ARE C DISTINCT AND IN ASCENDING ORDER IJ = -1 NMKM1 = NMK - 1 IREND = NMKM1 DO 70 I=1,NMKM1 E1 = ETA(I) + WK(I) E2 = ETA(I+1) + WK(I+1) IF (E1.LT.E2) GO TO 70 GO TO 90 70 CONTINUE IJ = 1 IREND = NMK C C TEST THAT THE SCHOENBERG-WHITNEY CONDITIONS ARE SATISFIED DO 80 I=1,NMK IPK = I + K TEMP = ETA(I) + WK(I) IF (TEMP.LE.X(I) .OR. TEMP.GE.X(IPK)) GO * TO 90 80 CONTINUE C THE SCHOENBERG-WHITNEY CONDITIONS ARE SATISFIED. C GO TO 160 C C REDUCE LAMBDA 90 LAMBDA = POINT5*LAMBDA IF (LAMBDA.LE.EPS) GO TO 250 C C IF DOING A CONTINUATION STEP THEN CONTINUE IF (ICONT.EQ.(-1)) GO TO 200 C C ELSE COMPUTE AN LAMBDA SUCH THAT IF IJ=-1 THE COMPONENTS OF C ETA +LAMBDA*(NEWTON CORRECTION) ARE DISTINCT AND IN C ASCENDING ORDER,OR IF IJ=1 THE COMPONENTS SATISFY THE C SCHOENBERG-WHITNEY CONDITIONS ITEND = 1 + IABS(INT(ALOG10(EPS)/ALOG10(TWO))) DO 130 IT=1,ITEND DO 110 IR=I,IREND IF (IJ.EQ.(-1)) GO TO 100 IRPK = IR + K TEMP = ETA(IR) + LAMBDA*WK(IR) IF (TEMP.LE.X(IR) .OR. * TEMP.GE.X(IRPK)) GO TO 120 GO TO 110 100 E1 = ETA(IR) + LAMBDA*WK(IR) E2 = ETA(IR+1) + LAMBDA*WK(IR+1) IF (E1.GE.E2) GO TO 120 110 CONTINUE GO TO 170 120 LAMBDA = POINT5*LAMBDA 130 CONTINUE C AT THIS STAGE LAMBDA IS .LE. EPS SO WE FINISH GO TO 250 C C FINISH ITERATION... C IF DOING A NEWTON STEP THEN RETURN 140 IF (ICONT.EQ.1) GO TO 260 C C ELSE IF ETA HAS NOT BEEN UPDATED THEN SET ERROR FLAG IF (IUP.EQ.(-1)) GO TO 250 C C OTHERWISE FINISH CONTINUATION AND START NEWTON STEP LAMBDA = ONE + LAMBDA ICONT = 1 C C STORE THE VECTOR ETA DO 150 I=1,NMK N3PI = N3 + I WK(N3PI) = ETA(I) 150 CONTINUE C END OF STORING ETA. GO TO 230 160 CONTINUE C C AT THIS STAGE NEWTONS METHOD HAS FAILED TO CONVERGE C AFTER IMAX ITERATIONS. GO TO 240 C C*************************************************************** C END OF NEWTON ITERATION. * C*************************************************************** C C*************************************************************** C CONTINUATION STEP * C*************************************************************** C C IF DOING A CONTINUATION STEP THEN CONTINUE 170 IF (ICONT.EQ.(-1)) GO TO 200 C C ELSE BEGIN CONTINUATION... C IF THIS IS THE FIRST CONTINUATION STEP THEN STORE THE C INITIAL APPROXIMATION IN WK(3*NMK+1),...,WK(4*NMK) IF (JCONT.GT.1) GO TO 190 DO 180 I=1,NMK N3PI = N3 + I WK(N3PI) = ETA(I) 180 CONTINUE C C STORE THE RESIDUALS, AT THE INITIAL APPROXIMATION, IN C WK(2*NMK+1),...,WK(3*NMK) 190 CALL RESID(IA, 1, NMK, N, X, NMK, WK(N3P1), CSTA, K, * WK(I1), WK(I2), N3, WK, N2, L, SS, ICONT, * LAMBDA) C C SET ICONT = -1 ICONT = -1 C C IF ETA HAS NOT BEEN UPDATED THEN BEGIN ANOTHER ITERATION 200 IF (IUP.EQ.(-1)) GO TO 230 C C ELSE,IF NECESSARY, RESET THE INITIAL APPROXIMATION IF (JCONT.EQ.1) GO TO 220 DO 210 I=1,NMK N3PI = N3 + I ETA(I) = WK(N3PI) 210 CONTINUE 220 IUP = -1 230 CONTINUE C C AT THIS STAGE THE CONTINUATION METHOD HAS FAILED TO C CONVERGE AFTER IMAX ITERATIONS. GO TO 240 C C*************************************************************** C END OF CONTINUATION METHOD. * C*************************************************************** C C*************************************************************** C SET THE ERROR RETURN FLAG: IFAIL * C*************************************************************** 240 IFAIL = 5 GO TO 260 250 IFAIL = 6 C END OF SETTING. C*************************************************************** 260 RETURN C*************************************************************** END SUBROUTINE RESID(IA, IARG, NR, N, X, NMK, ETA, CSTA, K, RES 10 * V, VINT, IW, WK, INC, L, SS, ICONT, LAMBDA) C ********************************************************* C PURPOSE C ********************************************************* C C THE PURPOSE OF THIS SUBROUTINE IS TWO-FOLD DEPENDING ON THE C VALUE OF IA. C C IF IA = +1 THEN THE SUBROUTINE COMPUTES: C A) THE INTEGERS L(J) WHERE X(L(J)).LE.ETA(J).LT.X(L(J)+1) C J=1,...,N-K. IT IS ASSUMED THAT THE KNOTS ETA(J), C J=1,...,N-K SATISFY THE SCHOENBERG-WHITNEY CONDITIONS. C C B) THE RESIDUALS AT ETA WHEN IARG = +1 C AND THE RESIDUALS AT ETA+(NEWTON CORRECTION) WHEN IARG = -1 C C C) THE SUM OF THE SQUARES OF THE RESIDUALS C WHERE THE RESIDUALS ARE: C F - LAMBDA*CSTA WHEN LAMBDA=ONE AND ICONT =+1 C AND C F - LAMBDA*CSTA - (ONE-LAMBDA)*CHAT WHEN ICONT =-1 C THE COMPONENTS OF CHAT ARE THE RESIDUALS AT THE INITIAL C APPROXIMATION. C C IN THIS CASE THE ARGUMENT NR IS THE NUMBER OF RESIDUALS. C THEREFORE IT SHOULD HAVE THE VALUE N-K. C C IF IA = -1 THEN THE SUBROUTINE COMPUTES: C C THE INTEGRAL OF THE QUANTITY MALPHA * PHIBAR WHERE MALPHA C IS THE AUGMENTED B-SPLINE AND PHIBAR IS A BANG-BANG C FUNCTION WHICH HAS BREAKPOINTS AT ETA(1),...,ETA(N-K), C AND WHICH HAS THE VALUE +1 ON THE INTERVAL C X(1).LE.X.LE.ETA(1). C C IN THIS CASE ONLY ONE INTEGRAL IS CALCULATED. THEREFORE C NR SHOULD HAVE THE VALUE +1. C C ********************************************************* REAL X(N), ETA(NMK), CSTA(NMK), V(K), VINT(K), WK(IW), * LAMBDA INTEGER L(NMK), P, Q, QP1 DATA ZERO, ONE, TWO /0.0,1.0,2.0/ SS = ZERO C ********************************************************* C IF IA=+1 THEN COMPUTE THE INTEGERS L(J),J=1,...,N-K C ********************************************************* IF (IA.EQ.(-1)) GO TO 40 DO 30 J=1,NMK IF (IARG.EQ.1) E1 = ETA(J) IF (IARG.EQ.(-1)) E1 = ETA(J) + WK(J) JPKM1 = J + K - 1 DO 10 JT=J,JPKM1 IF (E1.LT.X(JT+1)) GO TO 20 10 CONTINUE 20 L(J) = JT 30 CONTINUE C END OF COMPUTING L(J). C ********************************************************* C C ********************************************************* C BEGIN TO COMPUTE THE NR RESIDUALS C ********************************************************* 40 DO 330 I=1,NR IPK = I + K INCPI = INC + I KP1 = K + 1 IF (IA.EQ.1) GO TO 150 C C DETERMINE WHICH BREAKPOINTS ARE IN THE SUPPORT OF MALPHA C DO 50 Q=1,NMK IF (X(1).LE.ETA(Q)) GO TO 60 50 CONTINUE C C IF THIS COMMENT STATEMENT IS REACHED DURING EXECUTION THEN C THERE ARE NO BREAKPOINTS IN THE SUPPORT OF THE B-SPLINE MALPHA. C THEREFORE, THE INTEGRAL OF MALPHA * PHIBAR IS THE VALUE: C (-1)**IEXP*(1/K) WHERE IEXP = NMK GO TO 70 C 60 IF (ETA(Q).LT.X(KP1)) GO TO 80 C C LIKEWISE IF THIS STATEMENT IS REACHED THEN THERE ARE NO C BREAKPOINTS IN THE SUPPORT OF MALPHA AND THEREFORE THE INTEGRAL C OF MALPHA * PHIBAR IS THE VALUE: C (-1)**IEXP*(1/K) WHERE IEXP = Q - 1 70 WK(INCPI) = FLOAT((-1)**(IEXP))/FLOAT(K) RETURN C C THE SUPPORT OF THE B-SPLINE CONTAINS THE BREAKPOINT ETA(Q). C WE NOW DETERMINE THE REMAINING BREAKPOINTS IN THE SUPPORT. C 80 IF (Q.LT.NMK) GO TO 90 P = NMK GO TO 120 90 QP1 = Q + 1 DO 100 J=QP1,NMK IF (X(KP1).LE.ETA(J)) GO TO 110 100 CONTINUE J = NMK + 1 110 P = J - 1 C C THE SUPPORT OF MALPHA CONTAINS THE BREAKPOINTS ETA(Q),...,ETA(P). C C NOW COMPUTE THE KNOT INTERVALS WHERE THE DATA POINTS LIE 120 DO 140 J=Q,P E1 = ETA(J) IL = 1 IR = K + 2 130 MIDDLE = (IL+IR)/2 IF (E1.LT.X(MIDDLE)) IR = MIDDLE IF (E1.GE.X(MIDDLE)) IL = MIDDLE IF (IR.GT.IL+1) GO TO 130 L(J) = IL 140 CONTINUE C END OF COMPUTING THE KNOT INTERVALS. GO TO 200 C C IF IA=+1 THEN COMPUTE THE INTEGERS Q AND P C 150 J1 = MAX0(1,I-K+1) J2 = MIN0(I+1,NMK) C COMPUTE THE INTEGER Q SUCH THAT L(Q-1).LT.I.LE.L(Q) DO 160 Q=J1,NMK IF (I.LE.L(Q)) GO TO 170 160 CONTINUE C END OF COMPUTING Q. C C COMPUTE THE INTEGER P SUCH THAT L(P).LT.I+K.LE.L(P+1) 170 DO 180 J=J2,NMK IF (IPK.LE.L(J)) GO TO 190 180 CONTINUE J = NMK + 1 190 P = J - 1 C END OF COMPUTING P. C ********************************************************* C COMPUTE THE INTEGRALS AND STORE THE SUMS IN S1 AND S2 C ********************************************************* 200 DO 300 ISUM=1,2 SUM = ZERO IF (ISUM.EQ.1) GO TO 210 Q = Q + 1 IF (Q.GT.P) GO TO 290 210 DO 280 J=Q,P,2 IF (IARG.EQ.1) E1 = ETA(J) IF (IARG.EQ.(-1)) E1 = ETA(J) + WK(J) DO 220 IR=1,K V(IR) = ZERO VINT(IR) = ZERO 220 CONTINUE JINT = L(J) JJ1 = IPK - JINT E10 = E1 - X(JINT) E20 = X(JINT+1) - E1 V(1) = ONE/(X(JINT+1)-X(JINT)) VINT(1) = E10*V(1) NK = MIN0(K,JJ1) IF (NK.EQ.1) GO TO 240 DO 230 IJ=2,NK JINTPJ = JINT + IJ V(IJ) = (E10/(X(JINTPJ)-X(JINT)))* * V(IJ-1) VINT(IJ) = E10*V(IJ) 230 CONTINUE 240 IF (JINT.EQ.I) GO TO 270 MJ = JINT - I DO 260 IJ=1,MJ JINTMJ = JINT - IJ E30 = E1 - X(JINTMJ) V(1) = (E20/(X(JINT+1)-X(JINTMJ)))* * V(1) VINT(1) = VINT(1) + E30*V(1) NKJ = MIN0(K-IJ,JJ1) IF (NKJ.LE.1) GO TO 260 DO 250 LI=2,NKJ JINTPL = JINT + LI A1 = E30*V(LI-1) + * (X(JINTPL)-E1)*V(LI) A2 = X(JINTPL) - X(JINTMJ) V(LI) = A1/A2 VINT(LI) = VINT(LI) + E30*V(LI) 250 CONTINUE 260 CONTINUE 270 TERM = VINT(JJ1)/FLOAT(K) SUM = SUM + TERM 280 CONTINUE IF (ISUM.EQ.1) S1 = SUM 290 IF (ISUM.EQ.2) S2 = SUM 300 CONTINUE C END OF COMPUTING THE INTEGRALS. C ********************************************************* C C ********************************************************* C STORE THE RESIDUALS IN WK(INC+I),I=1,...,NR AND C COMPUTE THE SUM OF SQUARES IN SS C ********************************************************* SUM = FLOAT((-1)**(Q-2)) SUM = SUM*(S1-S2) IF (IA.EQ.1) GO TO 310 S1 = TWO*SUM + (FLOAT((-1)**P)/FLOAT(K)) GO TO 320 310 S1 = TWO*SUM + (FLOAT((-1)**P)/FLOAT(K)) - * LAMBDA*CSTA(I) N2PI = 2*NMK + I IF (ICONT.EQ.(-1)) S1 = S1 - (ONE-LAMBDA)*WK(N2PI) SS = SS + S1*S1 320 INCPI = INC + I WK(INCPI) = S1 330 CONTINUE C ********************************************************* C END OF COMPUTING THE RESIDUALS. C ********************************************************* RETURN END SUBROUTINE JAC(N, X, NMK, ETA, L, K, V, NK, AJ, Y, B) JAC 10 C ************************************************************* C PURPOSE C ************************************************************* C C GIVEN THE DATA POINTS: X(1)