1 DOUBLE PRECISION FUNCTION DBVALU(T,A,N,K,IDERIV,X,INBV,WORK) 2C***BEGIN PROLOGUE DBVALU 3C***DATE WRITTEN 800901 (YYMMDD) 4C***REVISION DATE 820801 (YYMMDD) 5C***REVISION HISTORY (YYMMDD) 6C 000330 Modified array declarations. (JEC) 7C 8C***CATEGORY NO. E3,K6 9C***KEYWORDS B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION, 10C SPLINE 11C***AUTHOR AMOS, D. E., (SNLA) 12C***PURPOSE Evaluates the B-representation of a B-spline at X for the 13C function value or any of its derivatives. 14C***DESCRIPTION 15C 16C Written by Carl de Boor and modified by D. E. Amos 17C 18C Reference 19C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. 20C 21C Abstract **** a double precision routine **** 22C DBVALU is the BVALUE function of the reference. 23C 24C DBVALU evaluates the B-representation (T,A,N,K) of a B-spline 25C at X for the function value on IDERIV=0 or any of its 26C derivatives on IDERIV=1,2,...,K-1. Right limiting values 27C (right derivatives) are returned except at the right end 28C point X=T(N+1) where left limiting values are computed. The 29C spline is defined on T(K) .LE. X .LE. T(N+1). DBVALU returns 30C a fatal error message when X is outside of this interval. 31C 32C To compute left derivatives or left limiting values at a 33C knot T(I), replace N by I-1 and set X=T(I), I=K+1,N+1. 34C 35C DBVALU calls DINTRV 36C 37C Description of Arguments 38C 39C Input T,A,X are double precision 40C T - knot vector of length N+K 41C A - B-spline coefficient vector of length N 42C N - number of B-spline coefficients 43C N = sum of knot multiplicities-K 44C K - order of the B-spline, K .GE. 1 45C IDERIV - order of the derivative, 0 .LE. IDERIV .LE. K-1 46C IDERIV = 0 returns the B-spline value 47C X - argument, T(K) .LE. X .LE. T(N+1) 48C INBV - an initialization parameter which must be set 49C to 1 the first time DBVALU is called. 50C 51C Output WORK,DBVALU are double precision 52C INBV - INBV contains information for efficient process- 53C ing after the initial call and INBV must not 54C be changed by the user. Distinct splines require 55C distinct INBV parameters. 56C WORK - work vector of length 3*K. 57C DBVALU - value of the IDERIV-th derivative at X 58C 59C Error Conditions 60C An improper input is a fatal error 61C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 62C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 63C JUNE 1977, PP. 441-472. 64C***ROUTINES CALLED DINTRV,XERROR 65C***END PROLOGUE DBVALU 66C 67C 68 INTEGER I,IDERIV,IDERP1,IHI,IHMKMJ,ILO,IMK,IMKPJ, INBV, IPJ, 69 1 IP1, IP1MJ, J, JJ, J1, J2, K, KMIDER, KMJ, KM1, KPK, MFLAG, N 70 DOUBLE PRECISION A, FKMJ, T, WORK, X 71 DIMENSION T(*), A(N), WORK(*) 72C***FIRST EXECUTABLE STATEMENT DBVALU 73 DBVALU = 0.0D0 74 IF(K.LT.1) GO TO 102 75 IF(N.LT.K) GO TO 101 76 IF(IDERIV.LT.0 .OR. IDERIV.GE.K) GO TO 110 77 KMIDER = K - IDERIV 78C 79C *** FIND *I* IN (K,N) SUCH THAT T(I) .LE. X .LT. T(I+1) 80C (OR, .LE. T(I+1) IF T(I) .LT. T(I+1) = T(N+1)). 81 KM1 = K - 1 82 CALL DINTRV(T, N+1, X, INBV, I, MFLAG) 83 IF (X.LT.T(K)) GO TO 120 84 IF (MFLAG.EQ.0) GO TO 20 85 IF (X.GT.T(I)) GO TO 130 86 10 IF (I.EQ.K) GO TO 140 87 I = I - 1 88 IF (X.EQ.T(I)) GO TO 10 89C 90C *** DIFFERENCE THE COEFFICIENTS *IDERIV* TIMES 91C WORK(I) = AJ(I), WORK(K+I) = DP(I), WORK(K+K+I) = DM(I), I=1.K 92C 93 20 IMK = I - K 94 DO 30 J=1,K 95 IMKPJ = IMK + J 96 WORK(J) = A(IMKPJ) 97 30 CONTINUE 98 IF (IDERIV.EQ.0) GO TO 60 99 DO 50 J=1,IDERIV 100 KMJ = K - J 101 FKMJ = DBLE(FLOAT(KMJ)) 102 DO 40 JJ=1,KMJ 103 IHI = I + JJ 104 IHMKMJ = IHI - KMJ 105 WORK(JJ) = (WORK(JJ+1)-WORK(JJ))/(T(IHI)-T(IHMKMJ))*FKMJ 106 40 CONTINUE 107 50 CONTINUE 108C 109C *** COMPUTE VALUE AT *X* IN (T(I),(T(I+1)) OF IDERIV-TH DERIVATIVE, 110C GIVEN ITS RELEVANT B-SPLINE COEFF. IN AJ(1),...,AJ(K-IDERIV). 111 60 IF (IDERIV.EQ.KM1) GO TO 100 112 IP1 = I + 1 113 KPK = K + K 114 J1 = K + 1 115 J2 = KPK + 1 116 DO 70 J=1,KMIDER 117 IPJ = I + J 118 WORK(J1) = T(IPJ) - X 119 IP1MJ = IP1 - J 120 WORK(J2) = X - T(IP1MJ) 121 J1 = J1 + 1 122 J2 = J2 + 1 123 70 CONTINUE 124 IDERP1 = IDERIV + 1 125 DO 90 J=IDERP1,KM1 126 KMJ = K - J 127 ILO = KMJ 128 DO 80 JJ=1,KMJ 129 WORK(JJ) = (WORK(JJ+1)*WORK(KPK+ILO)+WORK(JJ) 130 1 *WORK(K+JJ))/(WORK(KPK+ILO)+WORK(K+JJ)) 131 ILO = ILO - 1 132 80 CONTINUE 133 90 CONTINUE 134 100 DBVALU = WORK(1) 135 RETURN 136C 137C 138 101 CONTINUE 139! CALL XERROR( ' DBVALU, N DOES NOT SATISFY N.GE.K',35,2,1) 140 print *, ' DBVALU, N DOES NOT SATISFY N.GE.K' 141 RETURN 142 102 CONTINUE 143! CALL XERROR( ' DBVALU, K DOES NOT SATISFY K.GE.1',35,2,1) 144 print *, ' DBVALU, K DOES NOT SATISFY K.GE.1' 145 RETURN 146 110 CONTINUE 147! CALL XERROR( ' DBVALU, IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K', 148 print *, ' DBVALU, IDERIV DOES NOT SATISFY 0.LE.IDERIV.LT.K' 149 150 RETURN 151 120 CONTINUE 152! CALL XERROR( ' DBVALU, X IS N0T GREATER THAN OR EQUAL TO T(K)' 153 print *, ' DBVALU, X IS N0T GREATER THAN OR EQUAL TO T(K)' 154 RETURN 155 130 CONTINUE 156* CALL XERROR( ' DBVALU, X IS NOT LESS THAN OR EQUAL TO T(N+1)', 157* 1 47, 2, 1) 158 print *, ' DBVALU, X IS NOT LESS THAN OR EQUAL TO T(N+1)' 159 RETURN 160 140 CONTINUE 161* CALL XERROR( ' DBVALU, A LEFT LIMITING VALUE CANN0T BE OBTAINED A 162* 1T T(K)', 58, 2, 1) 163 print *,' DBVALU, A LEFT LIMITING VALUE CANT BE OBTAINED AT T(K)' 164 RETURN 165 END 166 167 SUBROUTINE DINTRV(XT,LXT,X,ILO,ILEFT,MFLAG) 168C***BEGIN PROLOGUE DINTRV 169C***DATE WRITTEN 800901 (YYMMDD) 170C***REVISION DATE 820801 (YYMMDD) 171C***CATEGORY NO. E3,K6 172C***KEYWORDS B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION, 173C SPLINE 174C***AUTHOR AMOS, D. E., (SNLA) 175C***PURPOSE Computes the largest integer ILEFT in 1.LE.ILEFT.LE.LXT 176C such that XT(ILEFT).LE.X where XT(*) is a subdivision of 177C the X interval. 178C***DESCRIPTION 179C 180C Written by Carl de Boor and modified by D. E. Amos 181C 182C Reference 183C SIAM J. Numerical Analysis, 14, No. 3, June 1977, pp.441-472. 184C 185C Abstract **** a double precision routine **** 186C DINTRV is the INTERV routine of the reference. 187C 188C DINTRV computes the largest integer ILEFT in 1 .LE. ILEFT .LE. 189C LXT such that XT(ILEFT) .LE. X where XT(*) is a subdivision of 190C the X interval. Precisely, 191C 192C X .LT. XT(1) 1 -1 193C if XT(I) .LE. X .LT. XT(I+1) then ILEFT=I , MFLAG=0 194C XT(LXT) .LE. X LXT 1, 195C 196C That is, when multiplicities are present in the break point 197C to the left of X, the largest index is taken for ILEFT. 198C 199C Description of Arguments 200C 201C Input XT,X are double precision 202C XT - XT is a knot or break point vector of length LXT 203C LXT - length of the XT vector 204C X - argument 205C ILO - an initialization parameter which must be set 206C to 1 the first time the spline array XT is 207C processed by DINTRV. 208C 209C Output 210C ILO - ILO contains information for efficient process- 211C ing after the initial call and ILO must not be 212C changed by the user. Distinct splines require 213C distinct ILO parameters. 214C ILEFT - largest integer satisfying XT(ILEFT) .LE. X 215C MFLAG - signals when X lies out of bounds 216C 217C Error Conditions 218C None 219C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 220C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 221C JUNE 1977, PP. 441-472. 222C***ROUTINES CALLED (NONE) 223C***END PROLOGUE DINTRV 224C 225C 226 INTEGER IHI, ILEFT, ILO, ISTEP, LXT, MFLAG, MIDDLE 227 DOUBLE PRECISION X, XT 228 DIMENSION XT(LXT) 229C***FIRST EXECUTABLE STATEMENT DINTRV 230 IHI = ILO + 1 231 IF (IHI.LT.LXT) GO TO 10 232 IF (X.GE.XT(LXT)) GO TO 110 233 IF (LXT.LE.1) GO TO 90 234 ILO = LXT - 1 235 IHI = LXT 236C 237 10 IF (X.GE.XT(IHI)) GO TO 40 238 IF (X.GE.XT(ILO)) GO TO 100 239C 240C *** NOW X .LT. XT(IHI) . FIND LOWER BOUND 241 ISTEP = 1 242 20 IHI = ILO 243 ILO = IHI - ISTEP 244 IF (ILO.LE.1) GO TO 30 245 IF (X.GE.XT(ILO)) GO TO 70 246 ISTEP = ISTEP*2 247 GO TO 20 248 30 ILO = 1 249 IF (X.LT.XT(1)) GO TO 90 250 GO TO 70 251C *** NOW X .GE. XT(ILO) . FIND UPPER BOUND 252 40 ISTEP = 1 253 50 ILO = IHI 254 IHI = ILO + ISTEP 255 IF (IHI.GE.LXT) GO TO 60 256 IF (X.LT.XT(IHI)) GO TO 70 257 ISTEP = ISTEP*2 258 GO TO 50 259 60 IF (X.GE.XT(LXT)) GO TO 110 260 IHI = LXT 261C 262C *** NOW XT(ILO) .LE. X .LT. XT(IHI) . NARROW THE INTERVAL 263 70 MIDDLE = (ILO+IHI)/2 264 IF (MIDDLE.EQ.ILO) GO TO 100 265C NOTE. IT IS ASSUMED THAT MIDDLE = ILO IN CASE IHI = ILO+1 266 IF (X.LT.XT(MIDDLE)) GO TO 80 267 ILO = MIDDLE 268 GO TO 70 269 80 IHI = MIDDLE 270 GO TO 70 271C *** SET OUTPUT AND RETURN 272 90 MFLAG = -1 273 ILEFT = 1 274 RETURN 275 100 MFLAG = 0 276 ILEFT = ILO 277 RETURN 278 110 MFLAG = 1 279 ILEFT = LXT 280 RETURN 281 END 282 283 SUBROUTINE DBKNOT(X,N,K,T) 284C***BEGIN PROLOGUE DBKNOT 285C***REFER TO DB2INK,DB3INK 286C***ROUTINES CALLED (NONE) 287C***REVISION HISTORY (YYMMDD) 288C 000330 Modified array declarations. (JEC) 289C 290C***END PROLOGUE DBKNOT 291C 292C -------------------------------------------------------------------- 293C DBKNOT CHOOSES A KNOT SEQUENCE FOR INTERPOLATION OF ORDER K AT THE 294C DATA POINTS X(I), I=1,..,N. THE N+K KNOTS ARE PLACED IN THE ARRAY 295C T. K KNOTS ARE PLACED AT EACH ENDPOINT AND NOT-A-KNOT END 296C CONDITIONS ARE USED. THE REMAINING KNOTS ARE PLACED AT DATA POINTS 297C IF N IS EVEN AND BETWEEN DATA POINTS IF N IS ODD. THE RIGHTMOST 298C KNOT IS SHIFTED SLIGHTLY TO THE RIGHT TO INSURE PROPER INTERPOLATION 299C AT X(N) (SEE PAGE 350 OF THE REFERENCE). 300C DOUBLE PRECISION VERSION OF BKNOT. 301C -------------------------------------------------------------------- 302C 303C ------------ 304C DECLARATIONS 305C ------------ 306C 307C PARAMETERS 308C 309 INTEGER 310 * N, K 311 DOUBLE PRECISION 312 * X(N), T(*) 313C 314C LOCAL VARIABLES 315C 316 INTEGER 317 * I, J, IPJ, NPJ, IP1 318 DOUBLE PRECISION 319 * RNOT 320C 321C 322C ---------------------------- 323C PUT K KNOTS AT EACH ENDPOINT 324C ---------------------------- 325C 326C (SHIFT RIGHT ENPOINTS SLIGHTLY -- SEE PG 350 OF REFERENCE) 327 RNOT = X(N) + 0.10D0*( X(N)-X(N-1) ) 328 DO 110 J=1,K 329 T(J) = X(1) 330 NPJ = N + J 331 T(NPJ) = RNOT 332 110 CONTINUE 333C 334C -------------------------- 335C DISTRIBUTE REMAINING KNOTS 336C -------------------------- 337C 338 IF (MOD(K,2) .EQ. 1) GO TO 150 339C 340C CASE OF EVEN K -- KNOTS AT DATA POINTS 341C 342 I = (K/2) - K 343 JSTRT = K+1 344 DO 120 J=JSTRT,N 345 IPJ = I + J 346 T(J) = X(IPJ) 347 120 CONTINUE 348 GO TO 200 349C 350C CASE OF ODD K -- KNOTS BETWEEN DATA POINTS 351C 352 150 CONTINUE 353 I = (K-1)/2 - K 354 IP1 = I + 1 355 JSTRT = K + 1 356 DO 160 J=JSTRT,N 357 IPJ = I + J 358 T(J) = 0.50D0*( X(IPJ) + X(IPJ+1) ) 359 160 CONTINUE 360 200 CONTINUE 361C 362 RETURN 363 END 364 365 SUBROUTINE DBTPCF(X,N,FCN,LDF,NF,T,K,BCOEF,WORK) 366C***BEGIN PROLOGUE DBTPCF 367C***REFER TO DB2INK,DB3INK 368C***ROUTINES CALLED DBINTK,DBNSLV 369C***REVISION HISTORY (YYMMDD) 370C 000330 Modified array declarations. (JEC) 371C 372C***END PROLOGUE DBTPCF 373C 374C ----------------------------------------------------------------- 375C DBTPCF COMPUTES B-SPLINE INTERPOLATION COEFFICIENTS FOR NF SETS 376C OF DATA STORED IN THE COLUMNS OF THE ARRAY FCN. THE B-SPLINE 377C COEFFICIENTS ARE STORED IN THE ROWS OF BCOEF HOWEVER. 378C EACH INTERPOLATION IS BASED ON THE N ABCISSA STORED IN THE 379C ARRAY X, AND THE N+K KNOTS STORED IN THE ARRAY T. THE ORDER 380C OF EACH INTERPOLATION IS K. THE WORK ARRAY MUST BE OF LENGTH 381C AT LEAST 2*K*(N+1). 382C DOUBLE PRECISION VERSION OF BTPCF. 383C ----------------------------------------------------------------- 384C 385C ------------ 386C DECLARATIONS 387C ------------ 388C 389C PARAMETERS 390C 391 INTEGER 392 * N, LDF, K 393 DOUBLE PRECISION 394 * X(N), FCN(LDF,NF), T(*), BCOEF(NF,N), WORK(*) 395C 396C LOCAL VARIABLES 397C 398 INTEGER 399 * I, J, K1, K2, IQ, IW 400C 401C --------------------------------------------- 402C CHECK FOR NULL INPUT AND PARTITION WORK ARRAY 403C --------------------------------------------- 404C 405C***FIRST EXECUTABLE STATEMENT 406 IF (NF .LE. 0) GO TO 500 407 K1 = K - 1 408 K2 = K1 + K 409 IQ = 1 + N 410 IW = IQ + K2*N+1 411C 412C ----------------------------- 413C COMPUTE B-SPLINE COEFFICIENTS 414C ----------------------------- 415C 416C 417C FIRST DATA SET 418C 419 CALL DBINTK(X,FCN,T,N,K,WORK,WORK(IQ),WORK(IW)) 420 DO 20 I=1,N 421 BCOEF(1,I) = WORK(I) 422 20 CONTINUE 423C 424C ALL REMAINING DATA SETS BY BACK-SUBSTITUTION 425C 426 IF (NF .EQ. 1) GO TO 500 427 DO 100 J=2,NF 428 DO 50 I=1,N 429 WORK(I) = FCN(I,J) 430 50 CONTINUE 431 CALL DBNSLV(WORK(IQ),K2,N,K1,K1,WORK) 432 DO 60 I=1,N 433 BCOEF(J,I) = WORK(I) 434 60 CONTINUE 435 100 CONTINUE 436C 437C ---- 438C EXIT 439C ---- 440C 441 500 CONTINUE 442 RETURN 443 END 444 445 SUBROUTINE DBINTK(X,Y,T,N,K,BCOEF,Q,WORK) 446C***BEGIN PROLOGUE DBINTK 447C***DATE WRITTEN 800901 (YYMMDD) 448C***REVISION DATE 820801 (YYMMDD) 449C***REVISION HISTORY (YYMMDD) 450C 000330 Modified array declarations. (JEC) 451C 452C***CATEGORY NO. E1A 453C***KEYWORDS B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION, 454C SPLINE 455C***AUTHOR AMOS, D. E., (SNLA) 456C***PURPOSE Produces the B-spline coefficients, BCOEF, of the 457C B-spline of order K with knots T(I), I=1,...,N+K, which 458C takes on the value Y(I) at X(I), I=1,...,N. 459C***DESCRIPTION 460C 461C Written by Carl de Boor and modified by D. E. Amos 462C 463C References 464C 465C A Practical Guide to Splines by C. de Boor, Applied 466C Mathematics Series 27, Springer, 1979. 467C 468C Abstract **** a double precision routine **** 469C 470C DBINTK is the SPLINT routine of the reference. 471C 472C DBINTK produces the B-spline coefficients, BCOEF, of the 473C B-spline of order K with knots T(I), I=1,...,N+K, which 474C takes on the value Y(I) at X(I), I=1,...,N. The spline or 475C any of its derivatives can be evaluated by calls to DBVALU. 476C 477C The I-th equation of the linear system A*BCOEF = B for the 478C coefficients of the interpolant enforces interpolation at 479C X(I), I=1,...,N. Hence, B(I) = Y(I), for all I, and A is 480C a band matrix with 2K-1 bands if A is invertible. The matrix 481C A is generated row by row and stored, diagonal by diagonal, 482C in the rows of Q, with the main diagonal going into row K. 483C The banded system is then solved by a call to DBNFAC (which 484C constructs the triangular factorization for A and stores it 485C again in Q), followed by a call to DBNSLV (which then 486C obtains the solution BCOEF by substitution). DBNFAC does no 487C pivoting, since the total positivity of the matrix A makes 488C this unnecessary. The linear system to be solved is 489C (theoretically) invertible if and only if 490C T(I) .LT. X(I) .LT. T(I+K), for all I. 491C Equality is permitted on the left for I=1 and on the right 492C for I=N when K knots are used at X(1) or X(N). Otherwise, 493C violation of this condition is certain to lead to an error. 494C 495C DBINTK calls DBSPVN, DBNFAC, DBNSLV, XERROR 496C 497C Description of Arguments 498C 499C Input X,Y,T are double precision 500C X - vector of length N containing data point abscissa 501C in strictly increasing order. 502C Y - corresponding vector of length N containing data 503C point ordinates. 504C T - knot vector of length N+K 505C Since T(1),..,T(K) .LE. X(1) and T(N+1),..,T(N+K) 506C .GE. X(N), this leaves only N-K knots (not nec- 507C essarily X(I) values) interior to (X(1),X(N)) 508C N - number of data points, N .GE. K 509C K - order of the spline, K .GE. 1 510C 511C Output BCOEF,Q,WORK are double precision 512C BCOEF - a vector of length N containing the B-spline 513C coefficients 514C Q - a work vector of length (2*K-1)*N, containing 515C the triangular factorization of the coefficient 516C matrix of the linear system being solved. The 517C coefficients for the interpolant of an 518C additional data set (X(I),yY(I)), I=1,...,N 519C with the same abscissa can be obtained by loading 520C YY into BCOEF and then executing 521C CALL DBNSLV(Q,2K-1,N,K-1,K-1,BCOEF) 522C WORK - work vector of length 2*K 523C 524C Error Conditions 525C Improper input is a fatal error 526C Singular system of equations is a fatal error 527C***REFERENCES C. DE BOOR, *A PRACTICAL GUIDE TO SPLINES*, APPLIED 528C MATHEMATICS SERIES 27, SPRINGER, 1979. 529C D.E. AMOS, *COMPUTATION WITH SPLINES AND B-SPLINES*, 530C SAND78-1968,SANDIA LABORATORIES,MARCH,1979. 531C***ROUTINES CALLED DBNFAC,DBNSLV,DBSPVN,XERROR 532C***END PROLOGUE DBINTK 533C 534C 535 INTEGER IFLAG, IWORK, K, N, I, ILP1MX, J, JJ, KM1, KPKM2, LEFT, 536 1 LENQ, NP1 537 DOUBLE PRECISION BCOEF(N), Y(N), Q(*), T(*), X(N), XI, WORK(*) 538C DIMENSION Q(2*K-1,N), T(N+K) 539C***FIRST EXECUTABLE STATEMENT DBINTK 540 IF(K.LT.1) GO TO 100 541 IF(N.LT.K) GO TO 105 542 JJ = N - 1 543 IF(JJ.EQ.0) GO TO 6 544 DO 5 I=1,JJ 545 IF(X(I).GE.X(I+1)) GO TO 110 546 5 CONTINUE 547 6 CONTINUE 548 NP1 = N + 1 549 KM1 = K - 1 550 KPKM2 = 2*KM1 551 LEFT = K 552C ZERO OUT ALL ENTRIES OF Q 553 LENQ = N*(K+KM1) 554 DO 10 I=1,LENQ 555 Q(I) = 0.0D0 556 10 CONTINUE 557C 558C *** LOOP OVER I TO CONSTRUCT THE N INTERPOLATION EQUATIONS 559 DO 50 I=1,N 560 XI = X(I) 561 ILP1MX = MIN0(I+K,NP1) 562C *** FIND LEFT IN THE CLOSED INTERVAL (I,I+K-1) SUCH THAT 563C T(LEFT) .LE. X(I) .LT. T(LEFT+1) 564C MATRIX IS SINGULAR IF THIS IS NOT POSSIBLE 565 LEFT = MAX0(LEFT,I) 566 IF (XI.LT.T(LEFT)) GO TO 80 567 20 IF (XI.LT.T(LEFT+1)) GO TO 30 568 LEFT = LEFT + 1 569 IF (LEFT.LT.ILP1MX) GO TO 20 570 LEFT = LEFT - 1 571 IF (XI.GT.T(LEFT+1)) GO TO 80 572C *** THE I-TH EQUATION ENFORCES INTERPOLATION AT XI, HENCE 573C A(I,J) = B(J,K,T)(XI), ALL J. ONLY THE K ENTRIES WITH J = 574C LEFT-K+1,...,LEFT ACTUALLY MIGHT BE NONZERO. THESE K NUMBERS 575C ARE RETURNED, IN BCOEF (USED FOR TEMP.STORAGE HERE), BY THE 576C FOLLOWING 577 30 CALL DBSPVN(T, K, K, 1, XI, LEFT, BCOEF, WORK, IWORK) 578C WE THEREFORE WANT BCOEF(J) = B(LEFT-K+J)(XI) TO GO INTO 579C A(I,LEFT-K+J), I.E., INTO Q(I-(LEFT+J)+2*K,(LEFT+J)-K) SINCE 580C A(I+J,J) IS TO GO INTO Q(I+K,J), ALL I,J, IF WE CONSIDER Q 581C AS A TWO-DIM. ARRAY , WITH 2*K-1 ROWS (SEE COMMENTS IN 582C DBNFAC). IN THE PRESENT PROGRAM, WE TREAT Q AS AN EQUIVALENT 583C ONE-DIMENSIONAL ARRAY (BECAUSE OF FORTRAN RESTRICTIONS ON 584C DIMENSION STATEMENTS) . WE THEREFORE WANT BCOEF(J) TO GO INTO 585C ENTRY 586C I -(LEFT+J) + 2*K + ((LEFT+J) - K-1)*(2*K-1) 587C = I-LEFT+1 + (LEFT -K)*(2*K-1) + (2*K-2)*J 588C OF Q . 589 JJ = I - LEFT + 1 + (LEFT-K)*(K+KM1) 590 DO 40 J=1,K 591 JJ = JJ + KPKM2 592 Q(JJ) = BCOEF(J) 593 40 CONTINUE 594 50 CONTINUE 595C 596C ***OBTAIN FACTORIZATION OF A , STORED AGAIN IN Q. 597 CALL DBNFAC(Q, K+KM1, N, KM1, KM1, IFLAG) 598 GO TO (60, 90), IFLAG 599C *** SOLVE A*BCOEF = Y BY BACKSUBSTITUTION 600 60 DO 70 I=1,N 601 BCOEF(I) = Y(I) 602 70 CONTINUE 603 CALL DBNSLV(Q, K+KM1, N, KM1, KM1, BCOEF) 604 RETURN 605C 606C 607 80 CONTINUE 608! CALL XERROR( ' DBINTK, SOME ABSCISSA WAS NOT IN THE SUPPORT OF TH 609! 1E CORRESPONDING BASIS FUNCTION AND THE SYSTEM IS SINGULAR.',109,2, 610! 21) 611 RETURN 612 90 CONTINUE 613! CALL XERROR( ' DBINTK, THE SYSTEM OF SOLVER DETECTS A SINGULAR SY 614! 1STEM ALTHOUGH THE THEORETICAL CONDITIONS FOR A SOLUTION WERE SATIS 615! 2FIED.',123,8,1) 616 RETURN 617 100 CONTINUE 618! CALL XERROR( ' DBINTK, K DOES NOT SATISFY K.GE.1', 35, 2, 1) 619 RETURN 620 105 CONTINUE 621! CALL XERROR( ' DBINTK, N DOES NOT SATISFY N.GE.K', 35, 2, 1) 622 RETURN 623 110 CONTINUE 624! CALL XERROR( ' DBINTK, X(I) DOES NOT SATISFY X(I).LT.X(I+1) FOR S 625! 1OME I', 57, 2, 1) 626 RETURN 627 END 628 629 SUBROUTINE DBNFAC(W,NROWW,NROW,NBANDL,NBANDU,IFLAG) 630C***BEGIN PROLOGUE DBNFAC 631C***REFER TO DBINT4,DBINTK 632C 633C DBNFAC is the BANFAC routine from 634C * A Practical Guide to Splines * by C. de Boor 635C 636C DBNFAC is a double precision routine 637C 638C Returns in W the LU-factorization (without pivoting) of the banded 639C matrix A of order NROW with (NBANDL + 1 + NBANDU) bands or diag- 640C onals in the work array W . 641C 642C ***** I N P U T ****** W is double precision 643C W.....Work array of size (NROWW,NROW) containing the interesting 644C part of a banded matrix A , with the diagonals or bands of A 645C stored in the rows of W , while columns of A correspond to 646C columns of W . This is the storage mode used in LINPACK and 647C results in efficient innermost loops. 648C Explicitly, A has NBANDL bands below the diagonal 649C + 1 (main) diagonal 650C + NBANDU bands above the diagonal 651C and thus, with MIDDLE = NBANDU + 1, 652C A(I+J,J) is in W(I+MIDDLE,J) for I=-NBANDU,...,NBANDL 653C J=1,...,NROW . 654C For example, the interesting entries of A (1,2)-banded matrix 655C of order 9 would appear in the first 1+1+2 = 4 rows of W 656C as follows. 657C 13 24 35 46 57 68 79 658C 12 23 34 45 56 67 78 89 659C 11 22 33 44 55 66 77 88 99 660C 21 32 43 54 65 76 87 98 661C 662C All other entries of W not identified in this way with an en- 663C try of A are never referenced . 664C NROWW.....Row dimension of the work array W . 665C must be .GE. NBANDL + 1 + NBANDU . 666C NBANDL.....Number of bands of A below the main diagonal 667C NBANDU.....Number of bands of A above the main diagonal . 668C 669C ***** O U T P U T ****** W is double precision 670C IFLAG.....Integer indicating success( = 1) or failure ( = 2) . 671C If IFLAG = 1, then 672C W.....contains the LU-factorization of A into a unit lower triangu- 673C lar matrix L and an upper triangular matrix U (both banded) 674C and stored in customary fashion over the corresponding entries 675C of A . This makes it possible to solve any particular linear 676C system A*X = B for X by a 677C CALL DBNSLV ( W, NROWW, NROW, NBANDL, NBANDU, B ) 678C with the solution X contained in B on return . 679C If IFLAG = 2, then 680C one of NROW-1, NBANDL,NBANDU failed to be nonnegative, or else 681C one of the potential pivots was found to be zero indicating 682C that A does not have an LU-factorization. This implies that 683C A is singular in case it is totally positive . 684C 685C ***** M E T H O D ****** 686C Gauss elimination W I T H O U T pivoting is used. The routine is 687C intended for use with matrices A which do not require row inter- 688C changes during factorization, especially for the T O T A L L Y 689C P O S I T I V E matrices which occur in spline calculations. 690C The routine should NOT be used for an arbitrary banded matrix. 691C***ROUTINES CALLED (NONE) 692C***END PROLOGUE DBNFAC 693C 694 INTEGER IFLAG, NBANDL, NBANDU, NROW, NROWW, I, IPK, J, JMAX, K, 695 1 KMAX, MIDDLE, MIDMK, NROWM1 696 DOUBLE PRECISION W(NROWW,NROW), FACTOR, PIVOT 697C 698C***FIRST EXECUTABLE STATEMENT DBNFAC 699 IFLAG = 1 700 MIDDLE = NBANDU + 1 701C W(MIDDLE,.) CONTAINS THE MAIN DIAGONAL OF A . 702 NROWM1 = NROW - 1 703 if (NROWM1 .lt. 0) then 704 goto 120 705 elseif (NROWM1 .eq. 0) then 706 goto 110 707 else 708 goto 10 709 endif 710 10 IF (NBANDL.GT.0) GO TO 30 711C A IS UPPER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO . 712 DO 20 I=1,NROWM1 713 IF (W(MIDDLE,I).EQ.0.0D0) GO TO 120 714 20 CONTINUE 715 GO TO 110 716 30 IF (NBANDU.GT.0) GO TO 60 717C A IS LOWER TRIANGULAR. CHECK THAT DIAGONAL IS NONZERO AND 718C DIVIDE EACH COLUMN BY ITS DIAGONAL . 719 DO 50 I=1,NROWM1 720 PIVOT = W(MIDDLE,I) 721 IF (PIVOT.EQ.0.0D0) GO TO 120 722 JMAX = MIN0(NBANDL,NROW-I) 723 DO 40 J=1,JMAX 724 W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT 725 40 CONTINUE 726 50 CONTINUE 727 RETURN 728C 729C A IS NOT JUST A TRIANGULAR MATRIX. CONSTRUCT LU FACTORIZATION 730 60 DO 100 I=1,NROWM1 731C W(MIDDLE,I) IS PIVOT FOR I-TH STEP . 732 PIVOT = W(MIDDLE,I) 733 IF (PIVOT.EQ.0.0D0) GO TO 120 734C JMAX IS THE NUMBER OF (NONZERO) ENTRIES IN COLUMN I 735C BELOW THE DIAGONAL . 736 JMAX = MIN0(NBANDL,NROW-I) 737C DIVIDE EACH ENTRY IN COLUMN I BELOW DIAGONAL BY PIVOT . 738 DO 70 J=1,JMAX 739 W(MIDDLE+J,I) = W(MIDDLE+J,I)/PIVOT 740 70 CONTINUE 741C KMAX IS THE NUMBER OF (NONZERO) ENTRIES IN ROW I TO 742C THE RIGHT OF THE DIAGONAL . 743 KMAX = MIN0(NBANDU,NROW-I) 744C SUBTRACT A(I,I+K)*(I-TH COLUMN) FROM (I+K)-TH COLUMN 745C (BELOW ROW I ) . 746 DO 90 K=1,KMAX 747 IPK = I + K 748 MIDMK = MIDDLE - K 749 FACTOR = W(MIDMK,IPK) 750 DO 80 J=1,JMAX 751 W(MIDMK+J,IPK) = W(MIDMK+J,IPK) - W(MIDDLE+J,I)*FACTOR 752 80 CONTINUE 753 90 CONTINUE 754 100 CONTINUE 755C CHECK THE LAST DIAGONAL ENTRY . 756 110 IF (W(MIDDLE,NROW).NE.0.0D0) RETURN 757 120 IFLAG = 2 758 RETURN 759 END 760 761 SUBROUTINE DBNSLV(W,NROWW,NROW,NBANDL,NBANDU,B) 762C***BEGIN PROLOGUE DBNSLV 763C***REFER TO DBINT4,DBINTK 764C 765C DBNSLV is the BANSLV routine from 766C * A Practical Guide to Splines * by C. de Boor 767C 768C DBNSLV is a double precision routine 769C 770C Companion routine to DBNFAC . It returns the solution X of the 771C linear system A*X = B in place of B , given the LU-factorization 772C for A in the work array W from DBNFAC. 773C 774C ***** I N P U T ****** W,B are DOUBLE PRECISION 775C W, NROWW,NROW,NBANDL,NBANDU.....Describe the LU-factorization of a 776C banded matrix A of order NROW as constructed in DBNFAC . 777C For details, see DBNFAC . 778C B.....Right side of the system to be solved . 779C 780C ***** O U T P U T ****** B is DOUBLE PRECISION 781C B.....Contains the solution X , of order NROW . 782C 783C ***** M E T H O D ****** 784C (With A = L*U, as stored in W,) the unit lower triangular system 785C L(U*X) = B is solved for Y = U*X, and Y stored in B . Then the 786C upper triangular system U*X = Y is solved for X . The calcul- 787C ations are so arranged that the innermost loops stay within columns. 788C***ROUTINES CALLED (NONE) 789C***END PROLOGUE DBNSLV 790C 791 INTEGER NBANDL, NBANDU, NROW, NROWW, I, J, JMAX, MIDDLE, NROWM1 792 DOUBLE PRECISION W(NROWW,NROW), B(NROW) 793C***FIRST EXECUTABLE STATEMENT DBNSLV 794 MIDDLE = NBANDU + 1 795 IF (NROW.EQ.1) GO TO 80 796 NROWM1 = NROW - 1 797 IF (NBANDL.EQ.0) GO TO 30 798C FORWARD PASS 799C FOR I=1,2,...,NROW-1, SUBTRACT RIGHT SIDE(I)*(I-TH COLUMN 800C OF L ) FROM RIGHT SIDE (BELOW I-TH ROW) . 801 DO 20 I=1,NROWM1 802 JMAX = MIN0(NBANDL,NROW-I) 803 DO 10 J=1,JMAX 804 B(I+J) = B(I+J) - B(I)*W(MIDDLE+J,I) 805 10 CONTINUE 806 20 CONTINUE 807C BACKWARD PASS 808C FOR I=NROW,NROW-1,...,1, DIVIDE RIGHT SIDE(I) BY I-TH DIAG- 809C ONAL ENTRY OF U, THEN SUBTRACT RIGHT SIDE(I)*(I-TH COLUMN 810C OF U) FROM RIGHT SIDE (ABOVE I-TH ROW). 811 30 IF (NBANDU.GT.0) GO TO 50 812C A IS LOWER TRIANGULAR . 813 DO 40 I=1,NROW 814 B(I) = B(I)/W(1,I) 815 40 CONTINUE 816 RETURN 817 50 I = NROW 818 60 B(I) = B(I)/W(MIDDLE,I) 819 JMAX = MIN0(NBANDU,I-1) 820 DO 70 J=1,JMAX 821 B(I-J) = B(I-J) - B(I)*W(MIDDLE-J,I) 822 70 CONTINUE 823 I = I - 1 824 IF (I.GT.1) GO TO 60 825 80 B(1) = B(1)/W(MIDDLE,1) 826 RETURN 827 END 828 829 SUBROUTINE DBSPVN(T,JHIGH,K,INDEX,X,ILEFT,VNIKX,WORK,IWORK) 830C***BEGIN PROLOGUE DBSPVN 831C***DATE WRITTEN 800901 (YYMMDD) 832C***REVISION DATE 820801 (YYMMDD) 833C***REVISION HISTORY (YYMMDD) 834C 000330 Modified array declarations. (JEC) 835C 836C***CATEGORY NO. E3,K6 837C***KEYWORDS B-SPLINE,DATA FITTING,DOUBLE PRECISION,INTERPOLATION, 838C SPLINE 839C***AUTHOR AMOS, D. E., (SNLA) 840C***PURPOSE Calculates the value of all (possibly) nonzero basis 841C functions at X. 842C***DESCRIPTION 843C 844C Written by Carl de Boor and modified by D. E. Amos 845C 846C Reference 847C SIAM J. Numerical Analysis, 14, No. 3, June, 1977, pp.441-472. 848C 849C Abstract **** a double precision routine **** 850C DBSPVN is the BSPLVN routine of the reference. 851C 852C DBSPVN calculates the value of all (possibly) nonzero basis 853C functions at X of order MAX(JHIGH,(J+1)*(INDEX-1)), where T(K) 854C .LE. X .LE. T(N+1) and J=IWORK is set inside the routine on 855C the first call when INDEX=1. ILEFT is such that T(ILEFT) .LE. 856C X .LT. T(ILEFT+1). A call to DINTRV(T,N+1,X,ILO,ILEFT,MFLAG) 857C produces the proper ILEFT. DBSPVN calculates using the basic 858C algorithm needed in DBSPVD. If only basis functions are 859C desired, setting JHIGH=K and INDEX=1 can be faster than 860C calling DBSPVD, but extra coding is required for derivatives 861C (INDEX=2) and DBSPVD is set up for this purpose. 862C 863C Left limiting values are set up as described in DBSPVD. 864C 865C Description of Arguments 866C 867C Input T,X are double precision 868C T - knot vector of length N+K, where 869C N = number of B-spline basis functions 870C N = sum of knot multiplicities-K 871C JHIGH - order of B-spline, 1 .LE. JHIGH .LE. K 872C K - highest possible order 873C INDEX - INDEX = 1 gives basis functions of order JHIGH 874C = 2 denotes previous entry with work, IWORK 875C values saved for subsequent calls to 876C DBSPVN. 877C X - argument of basis functions, 878C T(K) .LE. X .LE. T(N+1) 879C ILEFT - largest integer such that 880C T(ILEFT) .LE. X .LT. T(ILEFT+1) 881C 882C Output VNIKX, WORK are double precision 883C VNIKX - vector of length K for spline values. 884C WORK - a work vector of length 2*K 885C IWORK - a work parameter. Both WORK and IWORK contain 886C information necessary to continue for INDEX = 2. 887C When INDEX = 1 exclusively, these are scratch 888C variables and can be used for other purposes. 889C 890C Error Conditions 891C Improper input is a fatal error. 892C***REFERENCES C. DE BOOR, *PACKAGE FOR CALCULATING WITH B-SPLINES*, 893C SIAM JOURNAL ON NUMERICAL ANALYSIS, VOLUME 14, NO. 3, 894C JUNE 1977, PP. 441-472. 895C***ROUTINES CALLED XERROR 896C***END PROLOGUE DBSPVN 897C 898C 899 INTEGER ILEFT, IMJP1, INDEX, IPJ, IWORK, JHIGH, JP1, JP1ML, K, L 900 DOUBLE PRECISION T, VM, VMPREV, VNIKX, WORK, X 901C DIMENSION T(ILEFT+JHIGH) 902 DIMENSION T(*), VNIKX(K), WORK(*) 903C CONTENT OF J, DELTAM, DELTAP IS EXPECTED UNCHANGED BETWEEN CALLS. 904C WORK(I) = DELTAP(I), WORK(K+I) = DELTAM(I), I = 1,K 905C***FIRST EXECUTABLE STATEMENT DBSPVN 906 IF(K.LT.1) GO TO 90 907 IF(JHIGH.GT.K .OR. JHIGH.LT.1) GO TO 100 908 IF(INDEX.LT.1 .OR. INDEX.GT.2) GO TO 105 909 IF(X.LT.T(ILEFT) .OR. X.GT.T(ILEFT+1)) GO TO 110 910 GO TO (10, 20), INDEX 911 10 IWORK = 1 912 VNIKX(1) = 1.0D0 913 IF (IWORK.GE.JHIGH) GO TO 40 914C 915 20 IPJ = ILEFT + IWORK 916 WORK(IWORK) = T(IPJ) - X 917 IMJP1 = ILEFT - IWORK + 1 918 WORK(K+IWORK) = X - T(IMJP1) 919 VMPREV = 0.0D0 920 JP1 = IWORK + 1 921 DO 30 L=1,IWORK 922 JP1ML = JP1 - L 923 VM = VNIKX(L)/(WORK(L)+WORK(K+JP1ML)) 924 VNIKX(L) = VM*WORK(L) + VMPREV 925 VMPREV = VM*WORK(K+JP1ML) 926 30 CONTINUE 927 VNIKX(JP1) = VMPREV 928 IWORK = JP1 929 IF (IWORK.LT.JHIGH) GO TO 20 930C 931 40 RETURN 932C 933C 934 90 CONTINUE 935! CALL XERROR( ' DBSPVN, K DOES NOT SATISFY K.GE.1', 35, 2, 1) 936 RETURN 937 100 CONTINUE 938! CALL XERROR( ' DBSPVN, JHIGH DOES NOT SATISFY 1.LE.JHIGH.LE.K', 939! 1 48, 2, 1) 940 RETURN 941 105 CONTINUE 942! CALL XERROR( ' DBSPVN, INDEX IS NOT 1 OR 2',29,2,1) 943 RETURN 944 110 CONTINUE 945! CALL XERROR( ' DBSPVN, X DOES NOT SATISFY T(ILEFT).LE.X.LE.T(ILEF 946! 1T+1)', 56, 2, 1) 947 RETURN 948 END 949 950 DOUBLE PRECISION FUNCTION DB3VAL(XVAL,YVAL,ZVAL,IDX,IDY,IDZ, 951 * TX,TY,TZ,NX,NY,NZ,KX,KY,KZ,BCOEF,WORK) 952C***BEGIN PROLOGUE DB3VAL 953C***DATE WRITTEN 25 MAY 1982 954C***REVISION DATE 25 MAY 1982 955C***REVISION HISTORY (YYMMDD) 956C 000330 Modified array declarations. (JEC) 957C 958C***CATEGORY NO. E1A 959C***KEYWORDS INTERPOLATION, THREE-DIMENSIONS, GRIDDED DATA, SPLINES, 960C PIECEWISE POLYNOMIALS 961C***AUTHOR BOISVERT, RONALD, NBS 962C SCIENTIFIC COMPUTING DIVISION 963C NATIONAL BUREAU OF STANDARDS 964C WASHINGTON, DC 20234 965C***PURPOSE DB3VAL EVALUATES THE PIECEWISE POLYNOMIAL INTERPOLATING 966C FUNCTION CONSTRUCTED BY THE ROUTINE B3INK OR ONE OF ITS 967C PARTIAL DERIVATIVES. 968C DOUBLE PRECISION VERSION OF B3VAL. 969C***DESCRIPTION 970C 971C DB3VAL evaluates the tensor product piecewise polynomial 972C interpolant constructed by the routine DB3INK or one of its 973C derivatives at the point (XVAL,YVAL,ZVAL). To evaluate the 974C interpolant itself, set IDX=IDY=IDZ=0, to evaluate the first 975C partial with respect to x, set IDX=1,IDY=IDZ=0, and so on. 976C 977C DB3VAL returns 0.0D0 if (XVAL,YVAL,ZVAL) is out of range. That is, 978C XVAL.LT.TX(1) .OR. XVAL.GT.TX(NX+KX) .OR. 979C YVAL.LT.TY(1) .OR. YVAL.GT.TY(NY+KY) .OR. 980C ZVAL.LT.TZ(1) .OR. ZVAL.GT.TZ(NZ+KZ) 981C If the knots TX, TY, and TZ were chosen by DB3INK, then this is 982C equivalent to 983C XVAL.LT.X(1) .OR. XVAL.GT.X(NX)+EPSX .OR. 984C YVAL.LT.Y(1) .OR. YVAL.GT.Y(NY)+EPSY .OR. 985C ZVAL.LT.Z(1) .OR. ZVAL.GT.Z(NZ)+EPSZ 986C where EPSX = 0.1*(X(NX)-X(NX-1)), EPSY = 0.1*(Y(NY)-Y(NY-1)), and 987C EPSZ = 0.1*(Z(NZ)-Z(NZ-1)). 988C 989C The input quantities TX, TY, TZ, NX, NY, NZ, KX, KY, KZ, and BCOEF 990C should remain unchanged since the last call of DB3INK. 991C 992C 993C I N P U T 994C --------- 995C 996C XVAL Double precision scalar 997C X coordinate of evaluation point. 998C 999C YVAL Double precision scalar 1000C Y coordinate of evaluation point. 1001C 1002C ZVAL Double precision scalar 1003C Z coordinate of evaluation point. 1004C 1005C IDX Integer scalar 1006C X derivative of piecewise polynomial to evaluate. 1007C 1008C IDY Integer scalar 1009C Y derivative of piecewise polynomial to evaluate. 1010C 1011C IDZ Integer scalar 1012C Z derivative of piecewise polynomial to evaluate. 1013C 1014C TX Double precision 1D array (size NX+KX) 1015C Sequence of knots defining the piecewise polynomial in 1016C the x direction. (Same as in last call to DB3INK.) 1017C 1018C TY Double precision 1D array (size NY+KY) 1019C Sequence of knots defining the piecewise polynomial in 1020C the y direction. (Same as in last call to DB3INK.) 1021C 1022C TZ Double precision 1D array (size NZ+KZ) 1023C Sequence of knots defining the piecewise polynomial in 1024C the z direction. (Same as in last call to DB3INK.) 1025C 1026C NX Integer scalar 1027C The number of interpolation points in x. 1028C (Same as in last call to DB3INK.) 1029C 1030C NY Integer scalar 1031C The number of interpolation points in y. 1032C (Same as in last call to DB3INK.) 1033C 1034C NZ Integer scalar 1035C The number of interpolation points in z. 1036C (Same as in last call to DB3INK.) 1037C 1038C KX Integer scalar 1039C Order of polynomial pieces in x. 1040C (Same as in last call to DB3INK.) 1041C 1042C KY Integer scalar 1043C Order of polynomial pieces in y. 1044C (Same as in last call to DB3INK.) 1045C 1046C KZ Integer scalar 1047C Order of polynomial pieces in z. 1048C (Same as in last call to DB3INK.) 1049C 1050C BCOEF Double precision 2D array (size NX by NY by NZ) 1051C The B-spline coefficients computed by DB3INK. 1052C 1053C WORK Double precision 1D array (size KY*KZ+3*max(KX,KY,KZ)+KZ) 1054C A working storage array. 1055C 1056C***REFERENCES CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, 1057C SPRINGER-VERLAG, NEW YORK, 1978. 1058C***ROUTINES CALLED DINTRV,DBVALU 1059C***END PROLOGUE DB3VAL 1060C 1061C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 1062C 1063C MODIFICATION 1064C ------------ 1065C 1066C ADDED CHECK TO SEE IF X OR Y IS OUT OF RANGE, IF SO, RETURN 0.0 1067C 1068C R.F. BOISVERT, NIST 1069C 22 FEB 00 1070C 1071C<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> 1072C ------------ 1073C DECLARATIONS 1074C ------------ 1075C 1076C PARAMETERS 1077C 1078 INTEGER 1079 * IDX, IDY, IDZ, NX, NY, NZ, KX, KY, KZ 1080 DOUBLE PRECISION 1081 * XVAL, YVAL, ZVAL, TX(*), TY(*), TZ(*), BCOEF(NX,NY,NZ), 1082 * WORK(*) 1083C 1084C LOCAL VARIABLES 1085C 1086 INTEGER 1087 * ILOY, ILOZ, INBVX, INBV1, INBV2, LEFTY, LEFTZ, MFLAG, 1088 * KCOLY, KCOLZ, IZ, IZM1, IW, I, J, K 1089 DOUBLE PRECISION 1090 * DBVALU 1091C 1092 DATA ILOY /1/, ILOZ /1/, INBVX /1/ 1093C SAVE ILOY , ILOZ , INBVX 1094C 1095C 1096C***FIRST EXECUTABLE STATEMENT 1097 DB3VAL = 0.0D0 1098C NEXT STATEMENT - RFB MOD 1099 IF (XVAL.LT.TX(1) .OR. XVAL.GT.TX(NX+KX) .OR. 1100 + YVAL.LT.TY(1) .OR. YVAL.GT.TY(NY+KY) .OR. 1101 + ZVAL.LT.TZ(1) .OR. ZVAL.GT.TZ(NZ+KZ)) GO TO 100 1102 CALL DINTRV(TY,NY+KY,YVAL,ILOY,LEFTY,MFLAG) 1103 IF (MFLAG .NE. 0) GO TO 100 1104 CALL DINTRV(TZ,NZ+KZ,ZVAL,ILOZ,LEFTZ,MFLAG) 1105 IF (MFLAG .NE. 0) GO TO 100 1106 IZ = 1 + KY*KZ 1107 IW = IZ + KZ 1108 KCOLZ = LEFTZ - KZ 1109 I = 0 1110 DO 50 K=1,KZ 1111 KCOLZ = KCOLZ + 1 1112 KCOLY = LEFTY - KY 1113 DO 50 J=1,KY 1114 I = I + 1 1115 KCOLY = KCOLY + 1 1116 WORK(I) = DBVALU(TX,BCOEF(1,KCOLY,KCOLZ),NX,KX,IDX,XVAL, 1117 * INBVX,WORK(IW)) 1118 50 CONTINUE 1119 INBV1 = 1 1120 IZM1 = IZ - 1 1121 KCOLY = LEFTY - KY + 1 1122 DO 60 K=1,KZ 1123 I = (K-1)*KY + 1 1124 J = IZM1 + K 1125 WORK(J) = DBVALU(TY(KCOLY),WORK(I),KY,KY,IDY,YVAL, 1126 * INBV1,WORK(IW)) 1127 60 CONTINUE 1128 INBV2 = 1 1129 KCOLZ = LEFTZ - KZ + 1 1130 DB3VAL = DBVALU(TZ(KCOLZ),WORK(IZ),KZ,KZ,IDZ,ZVAL,INBV2, 1131 * WORK(IW)) 1132 100 CONTINUE 1133 RETURN 1134 END 1135 1136 SUBROUTINE DB3INK(X,NX,Y,NY,Z,NZ,FCN,LDF1,LDF2,KX,KY,KZ,TX,TY,TZ, 1137 * BCOEF,WORK,IFLAG) 1138C***BEGIN PROLOGUE DB3INK 1139C***DATE WRITTEN 25 MAY 1982 1140C***REVISION DATE 25 MAY 1982 1141C***REVISION HISTORY (YYMMDD) 1142C 000330 Modified array declarations. (JEC) 1143C 1144C***CATEGORY NO. E1A 1145C***KEYWORDS INTERPOLATION, THREE-DIMENSIONS, GRIDDED DATA, SPLINES, 1146C PIECEWISE POLYNOMIALS 1147C***AUTHOR BOISVERT, RONALD, NBS 1148C SCIENTIFIC COMPUTING DIVISION 1149C NATIONAL BUREAU OF STANDARDS 1150C WASHINGTON, DC 20234 1151C***PURPOSE DOUBLE PRECISION VERSION OF DB3INK 1152C DB3INK DETERMINES A PIECEWISE POLYNOMIAL FUNCTION THAT 1153C INTERPOLATES THREE-DIMENSIONAL GRIDDED DATA. USERS SPECIFY 1154C THE POLYNOMIAL ORDER (DEGREE+1) OF THE INTERPOLANT AND 1155C (OPTIONALLY) THE KNOT SEQUENCE. 1156C***DESCRIPTION 1157C 1158C DB3INK determines the parameters of a function that interpolates 1159C the three-dimensional gridded data (X(i),Y(j),Z(k),FCN(i,j,k)) for 1160C i=1,..,NX, j=1,..,NY, and k=1,..,NZ. The interpolating function and 1161C its derivatives may subsequently be evaluated by the function 1162C DB3VAL. 1163C 1164C The interpolating function is a piecewise polynomial function 1165C represented as a tensor product of one-dimensional B-splines. The 1166C form of this function is 1167C 1168C NX NY NZ 1169C S(x,y,z) = SUM SUM SUM a U (x) V (y) W (z) 1170C i=1 j=1 k=1 ij i j k 1171C 1172C where the functions U(i), V(j), and W(k) are one-dimensional B- 1173C spline basis functions. The coefficients a(i,j) are chosen so that 1174C 1175C S(X(i),Y(j),Z(k)) = FCN(i,j,k) for i=1,..,NX, j=1,..,NY, k=1,..,NZ 1176C 1177C Note that for fixed values of y and z S(x,y,z) is a piecewise 1178C polynomial function of x alone, for fixed values of x and z S(x,y, 1179C z) is a piecewise polynomial function of y alone, and for fixed 1180C values of x and y S(x,y,z) is a function of z alone. In one 1181C dimension a piecewise polynomial may be created by partitioning a 1182C given interval into subintervals and defining a distinct polynomial 1183C piece on each one. The points where adjacent subintervals meet are 1184C called knots. Each of the functions U(i), V(j), and W(k) above is a 1185C piecewise polynomial. 1186C 1187C Users of DB3INK choose the order (degree+1) of the polynomial 1188C pieces used to define the piecewise polynomial in each of the x, y, 1189C and z directions (KX, KY, and KZ). Users also may define their own 1190C knot sequence in x, y, and z separately (TX, TY, and TZ). If IFLAG= 1191C 0, however, DB3INK will choose sequences of knots that result in a 1192C piecewise polynomial interpolant with KX-2 continuous partial 1193C derivatives in x, KY-2 continuous partial derivatives in y, and KZ- 1194C 2 continuous partial derivatives in z. (KX knots are taken near 1195C each endpoint in x, not-a-knot end conditions are used, and the 1196C remaining knots are placed at data points if KX is even or at 1197C midpoints between data points if KX is odd. The y and z directions 1198C are treated similarly.) 1199C 1200C After a call to DB3INK, all information necessary to define the 1201C interpolating function are contained in the parameters NX, NY, NZ, 1202C KX, KY, KZ, TX, TY, TZ, and BCOEF. These quantities should not be 1203C altered until after the last call of the evaluation routine DB3VAL. 1204C 1205C 1206C I N P U T 1207C --------- 1208C 1209C X Double precision 1D array (size NX) 1210C Array of x abcissae. Must be strictly increasing. 1211C 1212C NX Integer scalar (.GE. 3) 1213C Number of x abcissae. 1214C 1215C Y Double precision 1D array (size NY) 1216C Array of y abcissae. Must be strictly increasing. 1217C 1218C NY Integer scalar (.GE. 3) 1219C Number of y abcissae. 1220C 1221C Z Double precision 1D array (size NZ) 1222C Array of z abcissae. Must be strictly increasing. 1223C 1224C NZ Integer scalar (.GE. 3) 1225C Number of z abcissae. 1226C 1227C FCN Double precision 3D array (size LDF1 by LDF2 by NY) 1228C Array of function values to interpolate. FCN(I,J,K) should 1229C contain the function value at the point (X(I),Y(J),Z(K)) 1230C 1231C LDF1 Integer scalar (.GE. NX) 1232C The actual first dimension of FCN used in the 1233C calling program. 1234C 1235C LDF2 Integer scalar (.GE. NY) 1236C The actual second dimension of FCN used in the calling 1237C program. 1238C 1239C KX Integer scalar (.GE. 2, .LT. NX) 1240C The order of spline pieces in x. 1241C (Order = polynomial degree + 1) 1242C 1243C KY Integer scalar (.GE. 2, .LT. NY) 1244C The order of spline pieces in y. 1245C (Order = polynomial degree + 1) 1246C 1247C KZ Integer scalar (.GE. 2, .LT. NZ) 1248C The order of spline pieces in z. 1249C (Order = polynomial degree + 1) 1250C 1251C 1252C I N P U T O R O U T P U T 1253C ----------------------------- 1254C 1255C TX Double precision 1D array (size NX+KX) 1256C The knots in the x direction for the spline interpolant. 1257C If IFLAG=0 these are chosen by DB3INK. 1258C If IFLAG=1 these are specified by the user. 1259C (Must be non-decreasing.) 1260C 1261C TY Double precision 1D array (size NY+KY) 1262C The knots in the y direction for the spline interpolant. 1263C If IFLAG=0 these are chosen by DB3INK. 1264C If IFLAG=1 these are specified by the user. 1265C (Must be non-decreasing.) 1266C 1267C TZ Double precision 1D array (size NZ+KZ) 1268C The knots in the z direction for the spline interpolant. 1269C If IFLAG=0 these are chosen by DB3INK. 1270C If IFLAG=1 these are specified by the user. 1271C (Must be non-decreasing.) 1272C 1273C 1274C O U T P U T 1275C ----------- 1276C 1277C BCOEF Double precision 3D array (size NX by NY by NZ) 1278C Array of coefficients of the B-spline interpolant. 1279C This may be the same array as FCN. 1280C 1281C 1282C M I S C E L L A N E O U S 1283C ------------------------- 1284C 1285C WORK Double precision 1D array (size NX*NY*NZ + max( 2*KX*(NX+1), 1286C 2*KY*(NY+1), 2*KZ*(NZ+1) ) 1287C Array of working storage. 1288C 1289C IFLAG Integer scalar. 1290C On input: 0 == knot sequence chosen by B2INK 1291C 1 == knot sequence chosen by user. 1292C On output: 1 == successful execution 1293C 2 == IFLAG out of range 1294C 3 == NX out of range 1295C 4 == KX out of range 1296C 5 == X not strictly increasing 1297C 6 == TX not non-decreasing 1298C 7 == NY out of range 1299C 8 == KY out of range 1300C 9 == Y not strictly increasing 1301C 10 == TY not non-decreasing 1302C 11 == NZ out of range 1303C 12 == KZ out of range 1304C 13 == Z not strictly increasing 1305C 14 == TY not non-decreasing 1306C 1307C***REFERENCES CARL DE BOOR, A PRACTICAL GUIDE TO SPLINES, 1308C SPRINGER-VERLAG, NEW YORK, 1978. 1309C CARL DE BOOR, EFFICIENT COMPUTER MANIPULATION OF TENSOR 1310C PRODUCTS, ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, 1311C VOL. 5 (1979), PP. 173-182. 1312C***ROUTINES CALLED DBTPCF,DBKNOT 1313C***END PROLOGUE DB3INK 1314C 1315C ------------ 1316C DECLARATIONS 1317C ------------ 1318C 1319C PARAMETERS 1320C 1321 INTEGER 1322 * NX, NY, NZ, LDF1, LDF2, KX, KY, KZ, IFLAG 1323 DOUBLE PRECISION 1324 * X(NX), Y(NY), Z(NZ), FCN(LDF1,LDF2,NZ), TX(*), TY(*), TZ(*), 1325 * BCOEF(NX,NY,NZ), WORK(*) 1326C 1327C LOCAL VARIABLES 1328C 1329 INTEGER 1330 * I, J, LOC, IW, NPK 1331C 1332C ----------------------- 1333C CHECK VALIDITY OF INPUT 1334C ----------------------- 1335C 1336C***FIRST EXECUTABLE STATEMENT 1337 IF ((IFLAG .LT. 0) .OR. (IFLAG .GT. 1)) GO TO 920 1338 IF (NX .LT. 3) GO TO 930 1339 IF (NY .LT. 3) GO TO 970 1340 IF (NZ .LT. 3) GO TO 1010 1341 IF ((KX .LT. 2) .OR. (KX .GE. NX)) GO TO 940 1342 IF ((KY .LT. 2) .OR. (KY .GE. NY)) GO TO 980 1343 IF ((KZ .LT. 2) .OR. (KZ .GE. NZ)) GO TO 1020 1344 DO 10 I=2,NX 1345 IF (X(I) .LE. X(I-1)) GO TO 950 1346 10 CONTINUE 1347 DO 20 I=2,NY 1348 IF (Y(I) .LE. Y(I-1)) GO TO 990 1349 20 CONTINUE 1350 DO 30 I=2,NZ 1351 IF (Z(I) .LE. Z(I-1)) GO TO 1030 1352 30 CONTINUE 1353 IF (IFLAG .EQ. 0) GO TO 70 1354 NPK = NX + KX 1355 DO 40 I=2,NPK 1356 IF (TX(I) .LT. TX(I-1)) GO TO 960 1357 40 CONTINUE 1358 NPK = NY + KY 1359 DO 50 I=2,NPK 1360 IF (TY(I) .LT. TY(I-1)) GO TO 1000 1361 50 CONTINUE 1362 NPK = NZ + KZ 1363 DO 60 I=2,NPK 1364 IF (TZ(I) .LT. TZ(I-1)) GO TO 1040 1365 60 CONTINUE 1366 70 CONTINUE 1367C 1368C ------------ 1369C CHOOSE KNOTS 1370C ------------ 1371C 1372 IF (IFLAG .NE. 0) GO TO 100 1373 CALL DBKNOT(X,NX,KX,TX) 1374 CALL DBKNOT(Y,NY,KY,TY) 1375 CALL DBKNOT(Z,NZ,KZ,TZ) 1376 100 CONTINUE 1377C 1378C ------------------------------- 1379C CONSTRUCT B-SPLINE COEFFICIENTS 1380C ------------------------------- 1381C 1382 IFLAG = 1 1383 IW = NX*NY*NZ + 1 1384C 1385C COPY FCN TO WORK IN PACKED FOR DBTPCF 1386 LOC = 0 1387 DO 510 K=1,NZ 1388 DO 510 J=1,NY 1389 DO 510 I=1,NX 1390 LOC = LOC + 1 1391 WORK(LOC) = FCN(I,J,K) 1392 510 CONTINUE 1393C 1394 CALL DBTPCF(X,NX,WORK,NX,NY*NZ,TX,KX,BCOEF,WORK(IW)) 1395 CALL DBTPCF(Y,NY,BCOEF,NY,NX*NZ,TY,KY,WORK,WORK(IW)) 1396 CALL DBTPCF(Z,NZ,WORK,NZ,NX*NY,TZ,KZ,BCOEF,WORK(IW)) 1397 GO TO 9999 1398C 1399C ----- 1400C EXITS 1401C ----- 1402C 1403 920 CONTINUE 1404! CALL XERRWV('DB3INK - IFLAG=I1 IS OUT OF RANGE.', 1405! * 35,2,1,1,IFLAG,I2,0,R1,R2) 1406 IFLAG = 2 1407 GO TO 9999 1408C 1409 930 CONTINUE 1410 IFLAG = 3 1411! CALL XERRWV('DB3INK - NX=I1 IS OUT OF RANGE.', 1412! * 32,IFLAG,1,1,NX,I2,0,R1,R2) 1413 GO TO 9999 1414C 1415 940 CONTINUE 1416 IFLAG = 4 1417! CALL XERRWV('DB3INK - KX=I1 IS OUT OF RANGE.', 1418! * 32,IFLAG,1,1,KX,I2,0,R1,R2) 1419 GO TO 9999 1420C 1421 950 CONTINUE 1422 IFLAG = 5 1423! CALL XERRWV('DB3INK - X ARRAY MUST BE STRICTLY INCREASING.', 1424! * 46,IFLAG,1,0,I1,I2,0,R1,R2) 1425 GO TO 9999 1426C 1427 960 CONTINUE 1428 IFLAG = 6 1429! CALL XERRWV('DB3INK - TX ARRAY MUST BE NON-DECREASING.', 1430! * 42,IFLAG,1,0,I1,I2,0,R1,R2) 1431 GO TO 9999 1432C 1433 970 CONTINUE 1434 IFLAG = 7 1435! CALL XERRWV('DB3INK - NY=I1 IS OUT OF RANGE.', 1436! * 32,IFLAG,1,1,NY,I2,0,R1,R2) 1437 GO TO 9999 1438C 1439 980 CONTINUE 1440 IFLAG = 8 1441! CALL XERRWV('DB3INK - KY=I1 IS OUT OF RANGE.', 1442! * 32,IFLAG,1,1,KY,I2,0,R1,R2) 1443 GO TO 9999 1444C 1445 990 CONTINUE 1446 IFLAG = 9 1447! CALL XERRWV('DB3INK - Y ARRAY MUST BE STRICTLY INCREASING.', 1448! * 46,IFLAG,1,0,I1,I2,0,R1,R2) 1449 GO TO 9999 1450C 1451 1000 CONTINUE 1452 IFLAG = 10 1453! CALL XERRWV('DB3INK - TY ARRAY MUST BE NON-DECREASING.', 1454! * 42,IFLAG,1,0,I1,I2,0,R1,R2) 1455 GO TO 9999 1456C 1457 1010 CONTINUE 1458 IFLAG = 11 1459! CALL XERRWV('DB3INK - NZ=I1 IS OUT OF RANGE.', 1460! * 32,IFLAG,1,1,NZ,I2,0,R1,R2) 1461 GO TO 9999 1462C 1463 1020 CONTINUE 1464 IFLAG = 12 1465! CALL XERRWV('DB3INK - KZ=I1 IS OUT OF RANGE.', 1466! * 32,IFLAG,1,1,KZ,I2,0,R1,R2) 1467 GO TO 9999 1468C 1469 1030 CONTINUE 1470 IFLAG = 13 1471! CALL XERRWV('DB3INK - Z ARRAY MUST BE STRICTLY INCREASING.', 1472! * 46,IFLAG,1,0,I1,I2,0,R1,R2) 1473 GO TO 9999 1474C 1475 1040 CONTINUE 1476 IFLAG = 14 1477! CALL XERRWV('DB3INK - TZ ARRAY MUST BE NON-DECREASING.', 1478! * 42,IFLAG,1,0,I1,I2,0,R1,R2) 1479 GO TO 9999 1480C 1481 9999 CONTINUE 1482 RETURN 1483 END 1484 1485