Changeset 5159 for LMDZ6/branches/Amaury_dev/libf/misc/lmdz_libmath_pch.f90
- Timestamp:
- Aug 2, 2024, 9:58:25 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/misc/lmdz_libmath_pch.f90
r5123 r5159 31 31 ! FTS 532-4275, (510) 422-4275 32 32 !***DESCRIPTION 33 ! 33 34 34 ! CHFEV: Cubic Hermite Function EValuator 35 ! 35 36 36 ! Evaluates the cubic polynomial determined by function values 37 37 ! F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points 38 38 ! XE(J), J=1(1)NE. 39 ! 39 40 40 ! ---------------------------------------------------------------------- 41 ! 41 42 42 ! Calling sequence: 43 ! 43 44 44 ! INTEGER NE, NEXT(2), IERR 45 45 ! REAL X1, X2, F1, F2, D1, D2, XE(NE), FE(NE) 46 ! 46 47 47 ! CALL CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR) 48 ! 48 49 49 ! Parameters: 50 ! 50 51 51 ! X1,X2 -- (input) endpoints of interval of definition of cubic. 52 52 ! (Error return if X1.EQ.X2 .) 53 ! 53 54 54 ! F1,F2 -- (input) values of function at X1 and X2, respectively. 55 ! 55 56 56 ! D1,D2 -- (input) values of derivative at X1 and X2, respectively. 57 ! 57 58 58 ! NE -- (input) number of evaluation points. (Error return if 59 59 ! NE.LT.1 .) 60 ! 60 61 61 ! XE -- (input) real array of points at which the function is to be 62 62 ! evaluated. If any of the XE are outside the interval 63 63 ! [X1,X2], a warning error is returned in NEXT. 64 ! 64 65 65 ! FE -- (output) real array of values of the cubic function defined 66 66 ! by X1,X2, F1,F2, D1,D2 at the points XE. 67 ! 67 68 68 ! NEXT -- (output) integer array indicating number of extrapolation 69 69 ! points: 70 70 ! NEXT(1) = number of evaluation points to left of interval. 71 71 ! NEXT(2) = number of evaluation points to right of interval. 72 ! 72 73 73 ! IERR -- (output) error flag. 74 74 ! Normal return: … … 78 78 ! IERR = -2 if X1.EQ.X2 . 79 79 ! (The FE-array has not been changed in either case.) 80 ! 80 81 81 !***REFERENCES (NONE) 82 82 !***ROUTINES CALLED XERMSG … … 92 92 !***END PROLOGUE CHFEV 93 93 ! Programming notes: 94 ! 94 95 95 ! To produce a double precision version, simply: 96 96 ! a. Change CHFEV to DCHFEV wherever it occurs, 97 97 ! b. Change the real declaration to double precision, and 98 98 ! c. Change the constant ZERO to double precision. 99 ! 99 100 100 ! DECLARE ARGUMENTS. 101 ! 101 102 102 INTEGER :: NE, NEXT(2), IERR 103 103 REAL :: X1, X2, F1, F2, D1, D2, XE(*), FE(*) 104 ! 104 105 105 ! DECLARE LOCAL VARIABLES. 106 ! 106 107 107 INTEGER :: I 108 108 REAL :: C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO 109 109 SAVE ZERO 110 110 DATA ZERO /0./ 111 ! 111 112 112 ! VALIDITY-CHECK ARGUMENTS. 113 ! 113 114 114 !***FIRST EXECUTABLE STATEMENT CHFEV 115 115 IF (NE < 1) GO TO 5001 116 116 H = X2 - X1 117 117 IF (H == ZERO) GO TO 5002 118 ! 118 119 119 ! INITIALIZE. 120 ! 120 121 121 IERR = 0 122 122 NEXT(1) = 0 … … 124 124 XMI = MIN(ZERO, H) 125 125 XMA = MAX(ZERO, H) 126 ! 126 127 127 ! COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1). 128 ! 128 129 129 DELTA = (F2 - F1) / H 130 130 DEL1 = (D1 - DELTA) / H … … 134 134 C3 = (DEL1 + DEL2) / H 135 135 ! (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.) 136 ! 136 137 137 ! EVALUATION LOOP. 138 ! 138 139 139 DO I = 1, NE 140 140 X = XE(I) - X1 … … 145 145 ! (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.) 146 146 END DO 147 ! 147 148 148 ! NORMAL RETURN. 149 ! 150 RETURN 151 ! 149 150 RETURN 151 152 152 ! ERROR RETURNS. 153 ! 153 154 154 5001 CONTINUE 155 155 ! NE.LT.1 RETURN. … … 158 158 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1) 159 159 RETURN 160 ! 160 161 161 5002 CONTINUE 162 162 ! X1.EQ.X2 RETURN. … … 184 184 ! FTS 532-4275, (510) 422-4275 185 185 !***DESCRIPTION 186 ! 186 187 187 ! PCHFE: Piecewise Cubic Hermite Function Evaluator 188 ! 188 189 189 ! Evaluates the cubic Hermite function defined by N, X, F, D at 190 190 ! the points XE(J), J=1(1)NE. 191 ! 191 192 192 ! To provide compatibility with PCHIM and PCHIC, includes an 193 193 ! increment between successive values of the F- and D-arrays. 194 ! 194 195 195 ! ---------------------------------------------------------------------- 196 ! 196 197 197 ! Calling sequence: 198 ! 198 199 199 ! PARAMETER (INCFD = ...) 200 200 ! INTEGER N, NE, IERR 201 201 ! REAL X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE) 202 202 ! LOGICAL SKIP 203 ! 203 204 204 ! CALL PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR) 205 ! 205 206 206 ! Parameters: 207 ! 207 208 208 ! N -- (input) number of data points. (Error return if N.LT.2 .) 209 ! 209 210 210 ! X -- (input) real array of independent variable values. The 211 211 ! elements of X must be strictly increasing: 212 212 ! X(I-1) .LT. X(I), I = 2(1)N. 213 213 ! (Error return if not.) 214 ! 214 215 215 ! F -- (input) real array of function values. F(1+(I-1)*INCFD) is 216 216 ! the value corresponding to X(I). 217 ! 217 218 218 ! D -- (input) real array of derivative values. D(1+(I-1)*INCFD) is 219 219 ! the value corresponding to X(I). 220 ! 220 221 221 ! INCFD -- (input) increment between successive values in F and D. 222 222 ! (Error return if INCFD.LT.1 .) 223 ! 223 224 224 ! SKIP -- (input/output) LOGICAL variable which should be set to 225 225 ! .TRUE. if the user wishes to skip checks for validity of … … 228 228 ! been performed (say, in PCHIM or PCHIC). 229 229 ! SKIP will be set to .TRUE. on normal return. 230 ! 230 231 231 ! NE -- (input) number of evaluation points. (Error return if 232 232 ! NE.LT.1 .) 233 ! 233 234 234 ! XE -- (input) real array of points at which the function is to be 235 235 ! evaluated. 236 ! 236 237 237 ! NOTES: 238 238 ! 1. The evaluation will be most efficient if the elements … … 243 243 ! values are extrapolated from the nearest extreme cubic, 244 244 ! and a warning error is returned. 245 ! 245 246 246 ! FE -- (output) real array of values of the cubic Hermite function 247 247 ! defined by N, X, F, D at the points XE. 248 ! 248 249 249 ! IERR -- (output) error flag. 250 250 ! Normal return: … … 261 261 ! NOTE: The above errors are checked in the order listed, 262 262 ! and following arguments have **NOT** been validated. 263 ! 263 264 264 !***REFERENCES (NONE) 265 265 !***ROUTINES CALLED CHFEV, XERMSG … … 275 275 !***END PROLOGUE PCHFE 276 276 ! Programming notes: 277 ! 277 278 278 ! 1. To produce a double precision version, simply: 279 279 ! a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they 280 280 ! occur, 281 281 ! b. Change the real declaration to double precision, 282 ! 282 283 283 ! 2. Most of the coding between the CALL to CHFEV and the end of 284 284 ! the IR-loop could be eliminated if it were permissible to 285 285 ! assume that XE is ordered relative to X. 286 ! 286 287 287 ! 3. CHFEV does not assume that X1 is less than X2. thus, it would 288 288 ! be possible to write a version of PCHFE that assumes a strict- 289 289 ! ly decreasing X-array by simply running the IR-loop backwards 290 290 ! (and reversing the order of appropriate tests). 291 ! 291 292 292 ! 4. The present code has a minor bug, which I have decided is not 293 293 ! worth the effort that would be required to fix it. … … 295 295 ! X(N-1), followed by points .GT.X(N), the extrapolation points 296 296 ! will be counted (at least) twice in the total returned in IERR. 297 ! 297 298 298 ! DECLARE ARGUMENTS. 299 ! 299 300 300 INTEGER :: N, INCFD, NE, IERR 301 301 REAL :: X(*), F(INCFD, *), D(INCFD, *), XE(*), FE(*) 302 302 LOGICAL :: SKIP 303 ! 303 304 304 ! DECLARE LOCAL VARIABLES. 305 ! 305 306 306 INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ 307 ! 307 308 308 ! VALIDITY-CHECK ARGUMENTS. 309 ! 309 310 310 !***FIRST EXECUTABLE STATEMENT PCHFE 311 311 IF (SKIP) GO TO 5 312 ! 312 313 313 IF (N<2) GO TO 5001 314 314 IF (INCFD<1) GO TO 5002 … … 316 316 IF (X(I)<=X(I - 1)) GO TO 5003 317 317 END DO 318 ! 318 319 319 ! FUNCTION DEFINITION IS OK, GO ON. 320 ! 320 321 321 5 CONTINUE 322 322 IF (NE<1) GO TO 5004 323 323 IERR = 0 324 324 SKIP = .TRUE. 325 ! 325 326 326 ! LOOP OVER INTERVALS. ( INTERVAL INDEX IS IL = IR-1 . ) 327 327 ! ( INTERVAL IS X(IL).LE.X.LT.X(IR) . ) … … 329 329 IR = 2 330 330 10 CONTINUE 331 ! 331 332 332 ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS. 333 ! 333 334 334 IF (JFIRST > NE) GO TO 5000 335 ! 335 336 336 ! LOCATE ALL POINTS IN INTERVAL. 337 ! 337 338 338 DO J = JFIRST, NE 339 339 IF (XE(J) >= X(IR)) GO TO 30 … … 341 341 J = NE + 1 342 342 GO TO 40 343 ! 343 344 344 ! HAVE LOCATED FIRST POINT BEYOND INTERVAL. 345 ! 345 346 346 30 CONTINUE 347 347 IF (IR == N) J = NE + 1 348 ! 348 349 349 40 CONTINUE 350 350 NJ = J - JFIRST 351 ! 351 352 352 ! SKIP EVALUATION IF NO POINTS IN INTERVAL. 353 ! 353 354 354 IF (NJ == 0) GO TO 50 355 ! 355 356 356 ! EVALUATE CUBIC AT XE(I), I = JFIRST (1) J-1 . 357 ! 357 358 358 ! ---------------------------------------------------------------- 359 359 CALL CHFEV (X(IR - 1), X(IR), F(1, IR - 1), F(1, IR), D(1, IR - 1), D(1, IR), & … … 361 361 ! ---------------------------------------------------------------- 362 362 IF (IERC < 0) GO TO 5005 363 ! 363 364 364 IF (NEXT(2) == 0) GO TO 42 365 365 ! IF (NEXT(2) .GT. 0) THEN 366 366 ! IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE 367 367 ! RIGHT OF X(IR). 368 ! 368 369 369 IF (IR < N) GO TO 41 370 370 ! IF (IR .EQ. N) THEN … … 379 379 ! ENDIF 380 380 42 CONTINUE 381 ! 381 382 382 IF (NEXT(1) == 0) GO TO 49 383 383 ! IF (NEXT(1) .GT. 0) THEN 384 384 ! IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE 385 385 ! LEFT OF X(IR-1). 386 ! 386 387 387 IF (IR > 2) GO TO 43 388 388 ! IF (IR .EQ. 2) THEN … … 394 394 ! XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST 395 395 ! EVALUATION INTERVAL. 396 ! 396 397 397 ! FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1). 398 398 DO I = JFIRST, J - 1 … … 402 402 ! IN CHFEV. 403 403 GO TO 5005 404 ! 404 405 405 45 CONTINUE 406 406 ! RESET J. (THIS WILL BE THE NEW JFIRST.) 407 407 J = I 408 ! 408 409 409 ! NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY. 410 410 DO I = 1, IR - 1 … … 412 412 END DO 413 413 ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1). 414 ! 414 415 415 47 CONTINUE 416 416 ! AT THIS POINT, EITHER XE(J) .LT. X(1) … … 422 422 ! ENDIF 423 423 49 CONTINUE 424 ! 424 425 425 JFIRST = J 426 ! 426 427 427 ! END OF IR-LOOP. 428 ! 428 429 429 50 CONTINUE 430 430 IR = IR + 1 431 431 IF (IR <= N) GO TO 10 432 ! 432 433 433 ! NORMAL RETURN. 434 ! 434 435 435 5000 CONTINUE 436 436 RETURN 437 ! 437 438 438 ! ERROR RETURNS. 439 ! 439 440 440 5001 CONTINUE 441 441 ! N.LT.2 RETURN. … … 444 444 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 445 445 RETURN 446 ! 446 447 447 5002 CONTINUE 448 448 ! INCFD.LT.1 RETURN. … … 451 451 1) 452 452 RETURN 453 ! 453 454 454 5003 CONTINUE 455 455 ! X-ARRAY NOT STRICTLY INCREASING. … … 458 458 , IERR, 1) 459 459 RETURN 460 ! 460 461 461 5004 CONTINUE 462 462 ! NE.LT.1 RETURN. … … 465 465 'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1) 466 466 RETURN 467 ! 467 468 468 5005 CONTINUE 469 469 ! ERROR RETURN FROM CHFEV. … … 662 662 !***AUTHOR Fritsch, F. N., (LLNL) 663 663 !***DESCRIPTION 664 ! 664 665 665 ! PCHDF: PCHIP Finite Difference Formula 666 ! 666 667 667 ! Uses a divided difference formulation to compute a K-point approx- 668 668 ! imation to the derivative at X(K) based on the data in X and S. 669 ! 669 670 670 ! Called by PCHCE and PCHSP to compute 3- and 4-point boundary 671 671 ! derivative approximations. 672 ! 672 673 673 ! ---------------------------------------------------------------------- 674 ! 674 675 675 ! On input: 676 676 ! K is the order of the desired derivative approximation. … … 682 682 ! S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1. 683 683 ! (Note that S need only be of length K-1.) 684 ! 684 685 685 ! On return: 686 686 ! S will be destroyed. … … 688 688 ! PCHDF will be set to the desired derivative approximation if 689 689 ! IERR=0 or to zero if IERR=-1. 690 ! 690 691 691 ! ---------------------------------------------------------------------- 692 ! 692 693 693 !***SEE ALSO PCHCE, PCHSP 694 694 !***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- … … 708 708 ! 930503 Improved purpose. (FNF) 709 709 !***END PROLOGUE PCHDF 710 ! 710 711 711 !**End 712 ! 712 713 713 ! DECLARE ARGUMENTS. 714 ! 714 715 715 INTEGER :: K, IERR 716 716 REAL :: X(K), S(K) 717 ! 717 718 718 ! DECLARE LOCAL VARIABLES. 719 ! 719 720 720 INTEGER :: I, J 721 721 REAL :: VALUE, ZERO 722 722 SAVE ZERO 723 723 DATA ZERO /0./ 724 ! 724 725 725 ! CHECK FOR LEGAL VALUE OF K. 726 ! 726 727 727 !***FIRST EXECUTABLE STATEMENT PCHDF 728 728 IF (K < 3) GO TO 5001 729 ! 729 730 730 ! COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL. 731 ! 731 732 732 DO J = 2, K - 1 733 733 DO I = 1, K - J … … 735 735 END DO 736 736 END DO 737 ! 737 738 738 ! EVALUATE DERIVATIVE AT X(K). 739 ! 739 740 740 VALUE = S(1) 741 741 DO I = 2, K - 1 742 742 VALUE = S(I) + VALUE * (X(K) - X(I)) 743 743 END DO 744 ! 744 745 745 ! NORMAL RETURN. 746 ! 746 747 747 IERR = 0 748 748 PCHDF = VALUE 749 749 RETURN 750 ! 750 751 751 ! ERROR RETURN. 752 ! 752 753 753 5001 CONTINUE 754 754 ! K.LT.3 RETURN. … … 776 776 ! FTS 532-4275, (510) 422-4275 777 777 !***DESCRIPTION 778 ! 778 779 779 ! PCHSP: Piecewise Cubic Hermite Spline 780 ! 780 781 781 ! Computes the Hermite representation of the cubic spline inter- 782 782 ! polant to the data given in X and F satisfying the boundary 783 783 ! conditions specified by IC and VC. 784 ! 784 785 785 ! To facilitate two-dimensional applications, includes an increment 786 786 ! between successive values of the F- and D-arrays. 787 ! 787 788 788 ! The resulting piecewise cubic Hermite function may be evaluated 789 789 ! by PCHFE or PCHFD. 790 ! 790 791 791 ! NOTE: This is a modified version of C. de Boor's cubic spline 792 792 ! routine CUBSPL. 793 ! 793 794 794 ! ---------------------------------------------------------------------- 795 ! 795 796 796 ! Calling sequence: 797 ! 797 798 798 ! PARAMETER (INCFD = ...) 799 799 ! INTEGER IC(2), N, NWK, IERR 800 800 ! REAL VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK) 801 ! 801 802 802 ! CALL PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR) 803 ! 803 804 804 ! Parameters: 805 ! 805 806 806 ! IC -- (input) integer array of length 2 specifying desired 807 807 ! boundary conditions: 808 808 ! IC(1) = IBEG, desired condition at beginning of data. 809 809 ! IC(2) = IEND, desired condition at end of data. 810 ! 810 811 811 ! IBEG = 0 to set D(1) so that the third derivative is con- 812 812 ! tinuous at X(2). This is the "not a knot" condition … … 823 823 ! 2. For the "natural" boundary condition, use IBEG=2 and 824 824 ! VC(1)=0. 825 ! 825 826 826 ! IEND may take on the same values as IBEG, but applied to 827 827 ! derivative at X(N). In case IEND = 1 or 2, the value is 828 828 ! given in VC(2). 829 ! 829 830 830 ! NOTES: 831 831 ! 1. An error return is taken if IEND is out of range. 832 832 ! 2. For the "natural" boundary condition, use IEND=2 and 833 833 ! VC(2)=0. 834 ! 834 835 835 ! VC -- (input) real array of length 2 specifying desired boundary 836 836 ! values, as indicated above. 837 837 ! VC(1) need be set only if IC(1) = 1 or 2 . 838 838 ! VC(2) need be set only if IC(2) = 1 or 2 . 839 ! 839 840 840 ! N -- (input) number of data points. (Error return if N.LT.2 .) 841 ! 841 842 842 ! X -- (input) real array of independent variable values. The 843 843 ! elements of X must be strictly increasing: 844 844 ! X(I-1) .LT. X(I), I = 2(1)N. 845 845 ! (Error return if not.) 846 ! 846 847 847 ! F -- (input) real array of dependent variable values to be inter- 848 848 ! polated. F(1+(I-1)*INCFD) is value corresponding to X(I). 849 ! 849 850 850 ! D -- (output) real array of derivative values at the data points. 851 851 ! These values will determine the cubic spline interpolant … … 854 854 ! D(1+(I-1)*INCFD), I=1(1)N. 855 855 ! No other entries in D are changed. 856 ! 856 857 857 ! INCFD -- (input) increment between successive values in F and D. 858 858 ! This argument is provided primarily for 2-D applications. 859 859 ! (Error return if INCFD.LT.1 .) 860 ! 860 861 861 ! WK -- (scratch) real array of working storage. 862 ! 862 863 863 ! NWK -- (input) length of work array. 864 864 ! (Error return if NWK.LT.2*N .) 865 ! 865 866 866 ! IERR -- (output) error flag. 867 867 ! Normal return: … … 882 882 ! (The D-array may have been changed in this case.) 883 883 ! ( Do **NOT** use it! ) 884 ! 884 885 885 !***REFERENCES Carl de Boor, A Practical Guide to Splines, Springer- 886 886 ! Verlag, New York, 1978, pp. 53-59. … … 899 899 !***END PROLOGUE PCHSP 900 900 ! Programming notes: 901 ! 901 902 902 ! To produce a double precision version, simply: 903 903 ! a. Change PCHSP to DPCHSP wherever it occurs, 904 904 ! b. Change the real declarations to double precision, and 905 905 ! c. Change the constants ZERO, HALF, ... to double precision. 906 ! 906 907 907 ! DECLARE ARGUMENTS. 908 ! 908 909 909 INTEGER :: IC(2), N, INCFD, NWK, IERR 910 910 REAL :: VC(2), X(*), F(INCFD, *), D(INCFD, *), WK(2, *) 911 ! 911 912 912 ! DECLARE LOCAL VARIABLES. 913 ! 913 914 914 INTEGER :: IBEG, IEND, INDEX, J, NM1 915 915 REAL :: G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO 916 916 SAVE ZERO, HALF, ONE, TWO, THREE 917 ! 917 918 918 DATA ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./ 919 ! 919 920 920 ! VALIDITY-CHECK ARGUMENTS. 921 ! 921 922 922 !***FIRST EXECUTABLE STATEMENT PCHSP 923 923 IF (N<2) GO TO 5001 … … 926 926 IF (X(J)<=X(J - 1)) GO TO 5003 927 927 END DO 928 ! 928 929 929 IBEG = IC(1) 930 930 IEND = IC(2) … … 933 933 IF ((IEND<0).OR.(IEND>4)) IERR = IERR - 2 934 934 IF (IERR<0) GO TO 5004 935 ! 935 936 936 ! FUNCTION DEFINITION IS OK -- GO ON. 937 ! 937 938 938 IF (NWK < 2 * N) GO TO 5007 939 ! 939 940 940 ! COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO, 941 941 ! COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.). … … 944 944 WK(2, J) = (F(1, J) - F(1, J - 1)) / WK(1, J) 945 945 END DO 946 ! 946 947 947 ! SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL. 948 ! 948 949 949 IF (IBEG>N) IBEG = 0 950 950 IF (IEND>N) IEND = 0 951 ! 951 952 952 ! SET UP FOR BOUNDARY CONDITIONS. 953 ! 953 954 954 IF ((IBEG==1).OR.(IBEG==2)) THEN 955 955 D(1, 1) = VC(1) … … 968 968 IBEG = 1 969 969 ENDIF 970 ! 970 971 971 IF ((IEND==1).OR.(IEND==2)) THEN 972 972 D(1, N) = VC(2) … … 985 985 IEND = 1 986 986 ENDIF 987 ! 987 988 988 ! --------------------( BEGIN CODING FROM CUBSPL )-------------------- 989 ! 989 990 990 ! **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF 991 991 ! F AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM- 992 992 ! INATION, WITH S(J) ENDING UP IN D(1,J), ALL J. 993 993 ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE. 994 ! 994 995 995 ! CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM 996 996 ! WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1) 997 ! 997 998 998 IF (IBEG == 0) THEN 999 999 IF (N == 2) THEN … … 1019 1019 D(1, 1) = THREE * WK(2, 2) - HALF * WK(1, 2) * D(1, 1) 1020 1020 ENDIF 1021 ! 1021 1022 1022 ! IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND 1023 1023 ! CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH 1024 1024 ! EQUATION READS WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J). 1025 ! 1025 1026 1026 NM1 = N - 1 1027 1027 IF (NM1 > 1) THEN … … 1034 1034 END DO 1035 1035 ENDIF 1036 ! 1036 1037 1037 ! CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM 1038 1038 ! (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N) 1039 ! 1039 1040 1040 ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK- 1041 1041 ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT 1042 1042 ! AT THIS POINT. 1043 1043 IF (IEND == 1) GO TO 30 1044 ! 1044 1045 1045 IF (IEND == 0) THEN 1046 1046 IF (N==2 .AND. IBEG==0) THEN … … 1073 1073 G = -ONE / WK(2, N - 1) 1074 1074 ENDIF 1075 ! 1075 1076 1076 ! COMPLETE FORWARD PASS OF GAUSS ELIMINATION. 1077 ! 1077 1078 1078 WK(2, N) = G * WK(1, N - 1) + WK(2, N) 1079 1079 IF (WK(2, N) == ZERO) GO TO 5008 1080 1080 D(1, N) = (G * D(1, N - 1) + D(1, N)) / WK(2, N) 1081 ! 1081 1082 1082 ! CARRY OUT BACK SUBSTITUTION 1083 ! 1083 1084 1084 30 CONTINUE 1085 1085 DO J = NM1, 1, -1 … … 1088 1088 END DO 1089 1089 ! --------------------( END CODING FROM CUBSPL )-------------------- 1090 ! 1090 1091 1091 ! NORMAL RETURN. 1092 ! 1093 RETURN 1094 ! 1092 1093 RETURN 1094 1095 1095 ! ERROR RETURNS. 1096 ! 1096 1097 1097 5001 CONTINUE 1098 1098 ! N.LT.2 RETURN. … … 1101 1101 'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1) 1102 1102 RETURN 1103 ! 1103 1104 1104 5002 CONTINUE 1105 1105 ! INCFD.LT.1 RETURN. … … 1108 1108 1) 1109 1109 RETURN 1110 ! 1110 1111 1111 5003 CONTINUE 1112 1112 ! X-ARRAY NOT STRICTLY INCREASING. … … 1115 1115 , IERR, 1) 1116 1116 RETURN 1117 ! 1117 1118 1118 5004 CONTINUE 1119 1119 ! IC OUT OF RANGE RETURN. … … 1121 1121 CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1) 1122 1122 RETURN 1123 ! 1123 1124 1124 5007 CONTINUE 1125 1125 ! NWK TOO SMALL RETURN. … … 1127 1127 CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1) 1128 1128 RETURN 1129 ! 1129 1130 1130 5008 CONTINUE 1131 1131 ! SINGULAR SYSTEM. … … 1136 1136 1) 1137 1137 RETURN 1138 ! 1138 1139 1139 5009 CONTINUE 1140 1140 ! ERROR RETURN FROM PCHDF.
Note: See TracChangeset
for help on using the changeset viewer.
