C ALGORITHM 770, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 23,NO. 2, June, 1997, P. 252--254. C #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Doc # Drivers # Src # This archive created: Mon Sep 1 10:20:59 1997 export PATH; PATH=/bin:$PATH if test ! -d 'Doc' then mkdir 'Doc' fi cd 'Doc' if test -f 'readme.doc' then echo shar: will not over-write existing file "'readme.doc'" else cat << \SHAR_EOF > 'readme.doc' **************************************** * * * BVSPIS package: some information * * * **************************************** The package is composed by the following 14 files: sdriver.for - ddriver.for source codes for drivers SBVSDR and DBVSDR illustrating, respectively, the use of SBVSPIS and DBVSPIS; sbvspis1.for - dbvspis1.for source codes for SBVSPIS and DBVSPIS single and double precision user callable routines; sbvspis2.for - dbvspis2.for source codes for SBVSPIS and DBVSPIS single and double precision internal routines; test1.inp first input file for SBVSPIS and DBVSPIS; test2.inp second input file for SBVSPIS and DBVSPIS; svermsf.out - dvermsf.out output files from SBVSDR and DBVSDR obtained running the code on a MS-DOS PC machine with Microsoft Fortran 77 compiler; svervax.out - dvervax.out output files from SBVSDR and DBVSDR obtained running the code on a VAX 6610 with DEC Fortran compiler; svercray.out - dvercray.out output files from SBVSDR and DBVSDR obtained running the code on a Cray C-90 with cf77 Cray compiler.  SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << \SHAR_EOF > 'data1' 10 .000000000E+00 .200000000E+00 .628318500E+00 .387785300E+00 .125663700E+01 .115105700E+01 .188495600E+01 .751056500E+00 .251327400E+01 .787785200E+00 .314159300E+01 -.200000100E+00 .376991100E+01 -.387785300E+00 .439823000E+01 -.115105700E+01 .502654800E+01 -.751056500E+00 .565486700E+01 -.787784900E+00 .628318500E+01 .200000200E+00 SHAR_EOF fi # end of overwriting check if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' DRESULTS.OUT THIS FILE CONTAINS THE COMPUTED RESULTS OF THE DRIVER PROGRAM DBVSDR, WHICH DEMONSTRATES THE USE OF THE SUBROUTINES DBVSIS, DBVSSC AND DBVSSE WITH SOME EXAMPLES EXAMPLE N. 1 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 -0.159087947 1 0.756827151 2 0.289082672 3 -0.289081974 4 -0.756825535 5 -0.935488313 6 -0.756826342 7 -0.289080577 8 0.289083319 9 0.756828209 10 2.387392074 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.000000000 0.200000000 1 0.628318500 0.387785300 2 1.256637000 1.151057000 3 1.884956000 0.751056500 4 2.513274000 0.787785200 5 3.141593000 -0.200000100 6 3.769911000 -0.387785300 7 4.398230000 -1.151057000 8 5.026548000 -0.751056500 9 5.654867000 -0.787784900 10 6.283185000 0.200000200 FURTHER EVALUATIONS VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 1.390000000 -0.472300291 -4.418754103 EXAMPLE N. 2 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 VALUES OF THE FIRST AND SECOND DERIVATIVES AT THE KNOTS I D(I) D2(I) 0 1.000000000 0.000000000 1 0.500000000 -0.250000000 2 0.000000000 -0.250000000 3 -0.541666667 -0.250000000 4 -0.833333333 0.000000000 5 -0.541666667 0.250000000 6 0.000000000 0.250000000 7 0.500000000 0.250000000 8 1.000000000 0.000000000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.000000000 0.000000000 0.000000000 1 0.800000000 0.642986667 -0.480000000 2 1.600000000 0.966048000 -0.586000000 3 2.400000000 0.970144000 -0.578000000 4 3.200000000 0.633102222 -0.613333333 5 4.000000000 0.000000000 0.000000000 6 4.800000000 -0.633102222 0.613333333 7 5.600000000 -0.970144000 0.578000000 8 6.400000000 -0.966048000 0.586000000 9 7.200000000 -0.642986667 0.480000000 10 8.000000000 0.000000000 0.000000000 EXAMPLE N. 3 VALUES OF THE ERROR FLAGS : ERRC=11 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 VALUES OF THE FIRST AND SECOND DERIVATIVES AT THE KNOTS I D(I) D2(I) 0 1.000000000 -4.390254282 1 0.000000000 0.411380511 2 0.000000000 -1.473300607 3 0.000000000 0.553123681 4 0.000000000 -0.697760587 5 0.000000000 0.000000000 6 0.000000000 0.000000000 7 0.000000000 1.473301230 8 0.000000000 -0.553124034 9 0.000000000 0.697757109 10 0.999998464 -4.999999693 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.000000000 0.200000000 1.000000000 1 0.418879000 0.387346154 0.062320677 2 0.837758000 0.548164656 1.776824196 3 1.256637000 1.151057000 0.000000000 4 1.675516000 0.835808140 -0.920285411 5 2.094395000 0.760405549 0.075775743 6 2.513274000 0.787785200 0.000000000 7 2.932153000 0.006557902 -2.318227415 8 3.351032000 -0.239411661 -0.442769193 9 3.769911000 -0.387785300 0.000000000 10 4.188790000 -0.984480728 -1.799682322 11 4.607669000 -1.061322842 0.934559864 12 5.026548000 -0.751056500 0.000000000 13 5.445427000 -0.777652602 -0.078018885 14 5.864306000 -0.624810338 1.802920032 15 6.283185000 0.200000200 0.999998464 FURTHER EVALUATIONS VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.000000000 -4.390254282 1 0.628318500 0.411380511 2 1.256637000 -1.473300607 3 1.884956000 0.553123681 4 2.513274000 -0.697760587 5 3.141593000 0.000000000 6 3.769911000 0.000000000 7 4.398230000 1.473301230 8 5.026548000 -0.553124034 9 5.654867000 0.697757109 10 6.283185000 -4.999999693 EXAMPLE N. 4 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 0.600000000 1 0.600000000 2 0.000000000 3 -0.600000000 4 -0.600000000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 -1.000000000 0.500000000 0.600000000 1 -0.866666667 0.580000000 0.600000000 2 -0.733333333 0.660000000 0.600000000 3 -0.600000000 0.740000000 0.600000000 4 -0.466666667 0.819970370 0.597333333 5 -0.333333333 0.896296296 0.533333333 6 -0.200000000 0.958400000 0.384000000 7 -0.066666667 0.994903704 0.149333333 8 0.066666667 0.994903704 -0.149333333 9 0.200000000 0.958400000 -0.384000000 10 0.333333333 0.896296296 -0.533333333 11 0.466666667 0.819970370 -0.597333333 12 0.600000000 0.740000000 -0.600000000 13 0.733333333 0.660000000 -0.600000000 14 0.866666667 0.580000000 -0.600000000 15 1.000000000 0.500000000 -0.600000000 EXAMPLE N. 5 VALUES OF THE INPUT CONSTRAINTS I CONSTR(I) 0 1 1 1 2 2 3 1 4 1 VALUES OF THE OUTPUT CONSTRAINTS I CONSTR(I) 0 1 1 1 2 2 3 1 4 1 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 1.000000000 1 1.000000000 2 0.500000000 3 -0.500000000 4 -1.000000000 5 -1.000000000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 1.270000000 1.291752449 1 -0.500000000 -0.500000000 2 5.700000000 -0.700000000 3 3.500000000 1.554687500 SHAR_EOF fi # end of overwriting check if test -f 'data2' then echo shar: will not over-write existing file "'data2'" else cat << \SHAR_EOF > 'data2' 8 0.0 0.0 1.0 0.75 2.0 1.0 3.0 0.75 4.0 0.0 5.0 -0.75 6.0 -1.0 7.0 -0.75 8.0 0.0 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' PROGRAM DBVSDR C ABSTRACT C C DBVSDR is a driver program which shows the use of the subroutines C DBVSIS, DBVSSC and DBVSSE with some examples. C C C DESCRIPTION C C DBVSDR consists of five examples which handle the execution of C DBVSSC (Boundary-Valued Shape-preserving Spline Computation) and C DBVSSE (Boundary-Valued Shape-preserving Spline Evaluation) or, when C convenient, the support routine DBVSIS, which simply collects the C calls to the previous ones, and cover the main features of the C package. The outputs are written to the external file RES. C For an easier reading of the program, all the output statements C are collected in the subroutine DWRTO. C The first and third examples make use of the data file data1 the C second of data2; obviously, they must be stored in the same C working directory as the executable version of DBVSDR. C The second and fourth examples make also use of external functions C which characterize the boundary conditions; although when the C non-separable boundary conditions are not imposed this function are C meaningless, vacuous function subprograms must be supplied for C linking purposes (for details see the description of the method and C of the input parameters BETA, BETAI, RHO and RHOI in DBVSIS). C C C DBVSDR makes use of: C C a) variables and arrays used as actual arguments for DBVSIS. C They have the same names of the corresponding dummy arguments C and their meaning is explained in the above subroutine; C C b) symbolic constant names, used for dimensioning the arrays. They C have the mnemonic form P-name, where P stands for PARAMETER and C 'name' is the name of a constant or variable which is used in C DBVSIS. So, their meaning can be easily obtained from the related C constant or variable 'name' described in the above subroutine; C C c) symbolic constants and variable names, used for input operations. C They will be explained at their first occurrence in executable C statements; C C d) function names used as actual arguments for DBVSIS. They have one C of the following forms DBn , DBnI , DRn , DRnI or DDUMFN. The C number n refers to the example in which the actual function is C used, B and R refer to functions 'beta' or 'rho' described in C DBVSIS and I , when present, to their inverse. DDUMFN is finally C used for linking purposes when 'beta' and/or 'rho' functions are C not referenced. C C The following are in class a) C C X,Y,NP,N,K,OPT,D0.DNP,D20,D2NP,CONSTR,EPS,KMAX,MAXSTP,XTAB,NTAB, C SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,Y0TAB,Y1TAB,Y2TAB,WORK,NWORK; C C in class b) C C PNP,PNTAB,PN,PCOMM,PPART,PNWORK; C C in class c) C C UN1,UN2,FNM1,FNM2; C C and in class d) C C DB2,DB2I,DR2,DR2I,DB4,DB4I,DDUMFN. C C DBVSDR calls the subroutines DBVSSC, DBVSSE and DWRTO C (at the entry points DWRT1A, DWRT1B, DWRT2, DWRT3A, DWRT3B, DWRT4, C DWRT5A, DWRT5B ). INTEGER PNP,PNTAB,PPART,PCOMM,PN,PNWORK PARAMETER (PNP=10,PNTAB=15,PPART=2,PCOMM=5,PN=6, * PNWORK=PCOMM+(PPART+7)*PNP+PN*(PN+11)/2+9) C The following statements define the symbolic names of constants used C for the input operations. INTEGER UN1,UN2 CHARACTER FNM1*9,FNM2*9 PARAMETER (UN1=1,UN2=2,FNM1='data1',FNM2='data2') INTEGER NP,N,K,OPT,CONSTR(0:PNP-1),KMAX,MAXSTP,NTAB,SBOPT, * Y0OPT,Y1OPT,Y2OPT,ERRE,ERRC,DIAGN(0:PNP-1),NWORK,I DOUBLE PRECISION X(0:PNP),Y(0:PNP),D0,DNP,D20,D2NP,EPS, * XTAB(0:PNTAB),D(0:PNP),D2(0:PNP), * Y0TAB(0:PNTAB),Y1TAB(0:PNTAB), * Y2TAB(0:PNTAB),WORK(1:PNWORK) DOUBLE PRECISION DB2,DB2I,DR2,DR2I,DB4,DB4I,DDUMFN EXTERNAL DB2,DB2I,DR2,DR2I,DB4,DB4I,DDUMFN C The following variables will be used in the output statements of C subroutine DWRTO. COMMON /DBVCOM/ NP,NTAB,ERRC,ERRE,X,Y,XTAB,Y0TAB,Y1TAB,Y2TAB, * D,D2,DIAGN,CONSTR C --------------------------------------------------------------------- C Call subroutine DWRTO for printing the heading in the output file. CALL DWRTO C --------------------------------------------------------------------- C Example n.1 . C C We want to construct a spline interpolating the set of data points C (x(i),y(i)) , i=0,1,...,10 , with x(i) = 2*i/20 and C y(i)= sin x(i) + ( cos 5x(i) )/5 , stored in the external file FNM1. C No constraints or boundary conditions are imposed on it and its C slopes are computed as the best approximation to the Bessel C estimates. The spline has continuity 1 and degree 3. Note that with C these selections we get the classical cubic Bessel interpolant. C We also want to estimate the spline at the set of tabulation points C xtab(i)=x(i) , i=0,...,10 . Then we want to evaluate the first and C the second derivative at the point 1.39 . C C C Read the data from the input file. It contains: C - the number NP of the data points, C - the coordinates (x(i),y(i)) of the i-th point, i=0,...,NP. OPEN (UNIT=UN1,FILE=FNM1,STATUS='OLD',ACCESS='SEQUENTIAL', * FORM='FORMATTED') READ (UNIT=UN1,FMT=*) NP READ (UNIT=UN1,FMT=*) (X(I),Y(I) , I=0,NP) CLOSE (UNIT=UN1,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=3 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: slopes as the best approximation to C the Bessel estimates, no boundary conditions, no constraints. OPT=210 C Set EPS, KMAX and MAXSTP to 0 for the automatic choice of their C values. EPS=0.D+00 KMAX=0 MAXSTP=0 C Call DBVSSC to compute the spline's parameters. CALL DBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,DDUMFN, * DDUMFN,DDUMFN,DDUMFN,KMAX,MAXSTP,ERRC,D,D2,DIAGN, * WORK,NWORK) C Set the options for evaluating the spline at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=10 DO 10 I=0,NTAB XTAB(I)=X(I) 10 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call DBVSSE for evaluating the spline. CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call DWRTO at the entry point DWRT1A for printing the results in the C output file. CALL DWRT1A C Set the options for evaluating the spline at the tabulation point. Y0OPT=0 Y1OPT=1 Y2OPT=1 C Define the new tabulation point. NTAB=0 XTAB(0)=0.139D+01 C Set SBOPT to 2 for selecting the binary search algorithm. SBOPT=2 C Call DBVSSE for evaluating the spline. CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call DWRTO at the entry point DWRT1B for printing the results in the C output file. CALL DWRT1B C --------------------------------------------------------------------- C Example n.2 . C C Given the set of points (x(i),y(i)) , i=0,...,8 where C {x(0),...,x(8)}={0.,1.,2.,3.,4.,5.,6.,7.,8.} and C {y(0),...,y(8)}={0.,0.75,1.,0.75,0.,-0.75,-1.,-0.75,0.} C stored in the external file FNM2, we want to construct a C co-monotone spline subject to the additional boundary condition C s'(8)=s'(0) , s''(8)=s''(0) . C We also ask the slopes to be the best approximation to the set of C third order monotonicity-preserving estimates described in DBVSSC C (see the input parameter OPT and the subroutine TODC). C The spline has degree 6 and continuity 2. C We evaluate the spline and its second derivative at the set of C tabulation points xtab(i), i=0,...,10 uniformly distributed in the C interval (x(0),x(8)). C We use in this case the support routine DBVSISP. C C Read the data from the input file. It contains: C - the number NP of the data points, C - the coordinates (x(i),y(i)) of the i-th point, i=0,...,8 . OPEN (UNIT=UN2,FILE=FNM2,STATUS='OLD', * ACCESS='SEQUENTIAL',FORM='FORMATTED') READ (UNIT=UN2,FMT=*) NP READ (UNIT=UN2,FMT=*) (X(I),Y(I),I=0,NP) CLOSE (UNIT=UN2,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=6 K=2 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: slopes as the best approximation to C the set of third order accurate estimates, non-separable boundary C conditions, monotonicity constraints. OPT=321 C Set EPS, KMAX and MAXSTP to 0 for the automatic choice of their C values. EPS=0.0D+00 KMAX=0 MAXSTP=0 C Set the options for evaluating the spline and its second derivative C at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=1 C Define the number and the set of tabulation points. NTAB=10 DO 20 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 20 CONTINUE C Set SBOPT to 1 for the sequential search algorithm. SBOPT=1 C Call DBVSIS for computing the spline. CALL DBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS, * DB2,DB2I,DR2,DR2I,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call DWRTO at the entry point DWRT2 for printing the results into C the output file. CALL DWRT2 C --------------------------------------------------------------------- C Example n.3 . C C Given the set of data points (x(i),y(i)) i=0,1,...,10 as in the C example 1, we want to construct a spline of degree 6 and continuity 2 C which is co-convex and subject to separable boundary conditions: C s'(x(0))=y'(x(0)) , s'(x(10))=y'(x(10)) , s''(x(0))=y''(x(0)) , C s''(x(10))=y''(x(10)) . The parameters are determined only by the C convexity constraints. We evaluate the spline and its first C derivative at the tabulation points xtab(i) , i=0,...,15 uniformly C distributed in the interval (x(0),x(10)) and then the second C derivative at the knots. C C C Read the data from the input file. OPEN (UNIT=UN1,FILE=FNM1,STATUS='OLD', * ACCESS='SEQUENTIAL',FORM='FORMATTED') READ (UNIT=UN1,FMT=*) NP READ (UNIT=UN1,FMT=*) (X(I),Y(I),I=0,NP) CLOSE (UNIT=UN1,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=6 K=2 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: the slopes are computed only from C the constraints, non-separable boundary conditions, C convexity constraints. OPT=132 C Set the boundary conditions. D0=COS(X(0))+SIN(5.0*X(0)) DNP=COS(X(NP))+SIN(5.0*X(NP)) D20=-SIN(X(0))-5.0*COS(5.0*X(0)) D2NP=-SIN(X(NP))-5.0*COS(5.0*X(NP)) C Set EPS, KMAX and MAXSTP for the automatic choice of their values. EPS=0.0D+00 KMAX=0 MAXSTP=0 C Call DBVSSC for computing the spline. CALL DBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,DDUMFN, * DDUMFN,DDUMFN,DDUMFN,KMAX,MAXSTP,ERRC,D,D2,DIAGN, * WORK,NWORK) C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=1 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=15 DO 30 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 30 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call DBVSSE for evaluating the spline. CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call DWRTO at the entry point DWRT3A for printing the results in the C output file. CALL DWRT3A C Set the options for evaluating the second derivative at the knots. Y0OPT=0 Y1OPT=0 Y2OPT=1 C Define the number and the set of tabulation points. NTAB=NP DO 40 I=0,NTAB XTAB(I)=X(I) 40 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call DBVSSE for evaluating the spline. CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call DWRTO at the entry point DWRT3B for printing the results in the C output file. CALL DWRT3B C --------------------------------------------------------------------- C Example n.4 . C C Given the 'Runge function': f(x) = 1/(1+x**2) , we want to construct C a co-monotone and co-convex spline which interpolates the function at C x(i) = -1+2*i/4 , i=0,1,...,4 . We also impose the additional C boundary conditions s'(1)=-s'(-1) and we compute the slopes as the C best approximation to the estimates s'(x(i)) = f'(x(i)) , i=0,...,4. C The spline has degree 3 and continuity 1. We evaluate the spline and C its first derivative at the tabulation points xtab(i), i=0,...,15 C uniformly distributed in the interval(x(0),x(4)). C It is used the support routine DBVSIS C C Compute the input data NP=4 DO 50 I=0,NP X(I)=-1.0+2.0*I/4.0 Y(I)=1.0/(1.0+X(I)**2) D(I)=-2.0*X(I)/((1.0+X(I)**2)**2) 50 CONTINUE C Set the degree and the class of continuity of the spline. N=3 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: best approximation to user supplied C derivatives, boundary conditions, monotonicity and convexity C constraints. OPT=423 C Set EPS, KMAX and MAXSTP EPS=1.0D-05 KMAX=5 MAXSTP=15 C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=1 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=15 DO 60 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 60 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call DBVSIS for computing the spline. CALL DBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,DB4,DB4I, * DDUMFN,DDUMFN,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call DWRTO at the entry point DWRT4 for printing the results into C the output file. CALL DWRT4 C --------------------------------------------------------------------- C Example n.5 . C C Given the points { (0,0),(1,1),(2,2),(3,2),(4,1),(5,0) } we want to C construct a spline of degree 4 and continuity 1 which is co-convex C in the interval (2,3) and co-monotone elsewhere. We also require the C separable boundary conditions s'(0)=1 , s'(5)=-1. The derivatives C are computed as the best approximation to Bessel estimates. C The spline is evaluated at four random points. C The routine DBVSIS is used. C C Compute the input data NP=5 X(0)=0.0D+00 X(1)=1.0D+00 X(2)=2.0D+00 X(3)=3.0D+00 X(4)=4.0D+00 X(5)=5.0D+00 Y(0)=0.0D+00 Y(1)=1.0D+00 Y(2)=2.0D+00 Y(3)=2.0D+00 Y(4)=1.0D+00 Y(5)=0.0D+00 C Set the degree and the class of continuity of the spline. N=4 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: Bessel like derivatives, C separable boundary conditions, user specified shape constraints. OPT=234 C Set the boundary conditions. D0=1.0D+00 DNP=-1.0D+00 C Select the shape constraints. CONSTR(0)=1 CONSTR(1)=1 CONSTR(2)=2 CONSTR(3)=1 CONSTR(4)=1 C Call DWRTO at the entry point DWRT5A for printing the input C constraints values. CALL DWRT5A C Set EPS, KMAX and MAXSTP for the automatic choice of their values. EPS=0. KMAX=0 MAXSTP=0 C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=3 XTAB(0)=0.127D+01 XTAB(1)=-0.5D+00 XTAB(2)=0.57D+01 XTAB(3)=0.35D+01 C Set SBOPT to 2 for selecting the binary search algorithm. SBOPT=2 C Call DBVSIS for computing the spline. CALL DBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,DDUMFN, * DDUMFN,DDUMFN,DDUMFN,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call DWRTO at the entry point DWRT5B for printing the results into C the output file. CALL DWRT5B STOP END C --------------------------------------------------------------------- SUBROUTINE DWRTO C DWRTO is an auxiliary routine which prints the results of the driver C program DBVSDR in the external file RES. All the parameters C are shared with the main program in the common area and their meaning C is explained in DBVSDR. C C Items of possible interest are: C C OUN : symbolic name of an integer constant, used for identifying C the logical output number; C FNAME : symbolic name of a character constant, used for identifying C the name of the output file. C C The subroutine is called at the entry points DWRTO, DWRT1A, DWRT1B, C DWRT2, DWRT3A, DWRT3B, DWRT4, DWRT5A, DWRT5B. INTEGER PNP,PNTAB,PN PARAMETER (PNP=10,PNTAB=15,PN=6) INTEGER OUN INTEGER I CHARACTER FNAME*12 PARAMETER (OUN=10,FNAME='RES') INTEGER NP,NTAB,CONSTR(0:PNP-1),ERRC,DIAGN(0:PNP-1),ERRE DOUBLE PRECISION X(0:PNP),Y(0:PNP),XTAB(0:PNTAB),Y0TAB(0:PNTAB), * Y1TAB(0:PNTAB),Y2TAB(0:PNTAB),D(0:PNP), * D2(0:PNP) COMMON /DBVCOM/ NP,NTAB,ERRC,ERRE,X,Y,XTAB,Y0TAB,Y1TAB,Y2TAB, * D,D2,DIAGN,CONSTR C Connect the output file FNAME to the logical unit number OUN. OPEN (UNIT=OUN,FILE=FNAME,STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') C Write the heading of the output file WRITE(UNIT=OUN,FMT=10) RETURN C Write the results of the example 1. ENTRY DWRT1A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 1' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=100) (I,XTAB(I),Y0TAB(I),I=0,NTAB) RETURN ENTRY DWRT1B WRITE (UNIT=OUN,FMT=25) 'FURTHER EVALUATIONS' WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=011) (I,XTAB(I),Y1TAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 2. ENTRY DWRT2 WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 2' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=45) WRITE (UNIT=OUN,FMT=55) (I,D(I),D2(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=101) (I,XTAB(I),Y0TAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 3. ENTRY DWRT3A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 3' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=45) WRITE (UNIT=OUN,FMT=55) (I,D(I),D2(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=110) (I,XTAB(I),Y0TAB(I),Y1TAB(I),I=0,NTAB) RETURN ENTRY DWRT3B WRITE (UNIT=OUN,FMT=25) 'FURTHER EVALUATIONS' WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=001) (I,XTAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 4. ENTRY DWRT4 WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 4' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=110) (I,XTAB(I),Y0TAB(I),Y1TAB(I),I=0,NTAB) RETURN C Write the results of the example 5. ENTRY DWRT5A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 5' WRITE (UNIT=OUN,FMT=70) WRITE (UNIT=OUN,FMT=90) (I,CONSTR(I),I=0,NP-1) RETURN ENTRY DWRT5B WRITE (UNIT=OUN,FMT=80) WRITE (UNIT=OUN,FMT=90) (I,CONSTR(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=100) (I,XTAB(I),Y0TAB(I),I=0,NTAB) C Disconnect the output file. CLOSE (UNIT=OUN,STATUS='KEEP') RETURN 10 FORMAT(//22X,'DRESULTS.OUT'// * ' THIS FILE CONTAINS THE COMPUTED RESULTS OF THE DRIVER'/ * ' PROGRAM DBVSDR, WHICH DEMONSTRATES THE USE OF THE'/ * ' SUBROUTINES DBVSIS, DBVSSC AND DBVSSE WITH SOME EXAMPLES') 20 FORMAT(//22X,A13/) 25 FORMAT(/22X,A19/) 30 FORMAT(/' VALUES OF THE ERROR FLAGS : ERRC=', * I2,' ; ERRE=',I2) 33 FORMAT(/' DIAGNOSTIC INFORMATION'//, * 10X,'I',7X,'DIAGN(I)'/) 36 FORMAT(9X,I2,11X,I1) 40 FORMAT(/' VALUES OF THE FIRST DERIVATIVES AT THE KNOTS'//, * 10X,'I',11X,'D(I)'/) 45 FORMAT(/' VALUES OF THE FIRST AND SECOND DERIVATIVES AT ', * 'THE KNOTS'//,10X,'I',11X,'D(I)',14X,'D2(I)'/) 50 FORMAT(9X,I2,5X,F13.9) 55 FORMAT(9X,I2,5X,F13.9,5X,F13.9) 60 FORMAT(/ * ' VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION', * ' POINTS'//, * 7X,'I',8X,'XTAB(I)',8X,'Y0TAB(I)',8X,'Y1TAB(I)',8X,'Y2TAB(I)'/) 70 FORMAT(/' VALUES OF THE INPUT CONSTRAINTS'//, * 10X,'I',7X,'CONSTR(I)'/) 80 FORMAT(/' VALUES OF THE OUTPUT CONSTRAINTS'//, * 10X,'I',7X,'CONSTR(I)'/) 90 FORMAT(9X,I2,11X,I1) 001 FORMAT(5X,I3,3X,F13.9,35X,F13.9) 011 FORMAT(5X,I3,3X,F13.9,19X,F13.9,3X,F13.9) 100 FORMAT(5X,I3,3X,F13.9,3X,F13.9) 101 FORMAT(5X,I3,3X,F13.9,3X,F13.9,19X,F13.9) 110 FORMAT(5X,I3,3X,F13.9,3X,F13.9,3X,F13.9) END C --------------------------------------------------------------------- C The following functions define the non-separable boundary conditions. C For linking purposes, these must be given even when the conditions C are not required. C In this case the 'vacuous function' DDUMFN is used. DOUBLE PRECISION FUNCTION DB2(X) DOUBLE PRECISION X DB2=X RETURN END DOUBLE PRECISION FUNCTION DB2I(Y) DOUBLE PRECISION Y DB2I=Y RETURN END DOUBLE PRECISION FUNCTION DR2(X) DOUBLE PRECISION X DR2=X RETURN END DOUBLE PRECISION FUNCTION DR2I(Y) DOUBLE PRECISION Y DR2I=Y RETURN END DOUBLE PRECISION FUNCTION DB4(X) DOUBLE PRECISION X DB4=-X RETURN END DOUBLE PRECISION FUNCTION DB4I(Y) DOUBLE PRECISION Y DB4I=-Y RETURN END DOUBLE PRECISION FUNCTION DDUMFN(T) DOUBLE PRECISION T DDUMFN=0.0D+00*T RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' SRESULTS.OUT THIS FILE CONTAINS THE COMPUTED RESULTS OF THE DRIVER PROGRAM SBVSDR, WHICH DEMONSTRATES THE USE OF THE SUBROUTINES SBVSIS, SBVSSC AND SBVSSE WITH SOME EXAMPLES EXAMPLE N. 1 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 -0.15909 1 0.75683 2 0.28908 3 -0.28908 4 -0.75683 5 -0.93549 6 -0.75683 7 -0.28908 8 0.28908 9 0.75683 10 2.38739 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.00000 0.20000 1 0.62832 0.38779 2 1.25664 1.15106 3 1.88496 0.75106 4 2.51327 0.78779 5 3.14159 -0.20000 6 3.76991 -0.38779 7 4.39823 -1.15106 8 5.02655 -0.75106 9 5.65487 -0.78778 10 6.28319 0.20000 FURTHER EVALUATIONS VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 1.39000 -0.47230 -4.41875 EXAMPLE N. 2 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 VALUES OF THE FIRST AND SECOND DERIVATIVES AT THE KNOTS I D(I) D2(I) 0 1.00000 0.00000 1 0.50000 -0.25000 2 0.00000 -0.25000 3 -0.54167 -0.25000 4 -0.83333 0.00000 5 -0.54167 0.25000 6 0.00000 0.25000 7 0.50000 0.25000 8 1.00000 0.00000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.00000 0.00000 0.00000 1 0.80000 0.64299 -0.48000 2 1.60000 0.96605 -0.58600 3 2.40000 0.97014 -0.57800 4 3.20000 0.63310 -0.61333 5 4.00000 0.00000 0.00000 6 4.80000 -0.63310 0.61333 7 5.60000 -0.97014 0.57800 8 6.40000 -0.96605 0.58600 9 7.20000 -0.64299 0.48000 10 8.00000 0.00000 0.00000 EXAMPLE N. 3 VALUES OF THE ERROR FLAGS : ERRC=11 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 9 0 VALUES OF THE FIRST AND SECOND DERIVATIVES AT THE KNOTS I D(I) D2(I) 0 1.00000 -4.39025 1 0.00000 0.41138 2 0.00000 -1.47330 3 0.00000 0.55312 4 0.00000 -0.69776 5 0.00000 0.00000 6 0.00000 0.00000 7 0.00000 1.47330 8 0.00000 -0.55312 9 0.00000 0.69776 10 1.00000 -5.00000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.00000 0.20000 1.00000 1 0.41888 0.38735 0.06232 2 0.83776 0.54816 1.77682 3 1.25664 1.15106 0.00000 4 1.67552 0.83581 -0.92029 5 2.09439 0.76041 0.07578 6 2.51327 0.78779 0.00000 7 2.93215 0.00656 -2.31823 8 3.35103 -0.23941 -0.44277 9 3.76991 -0.38779 0.00000 10 4.18879 -0.98448 -1.79968 11 4.60767 -1.06132 0.93456 12 5.02655 -0.75106 0.00000 13 5.44543 -0.77765 -0.07802 14 5.86431 -0.62481 1.80292 15 6.28319 0.20000 1.00000 FURTHER EVALUATIONS VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 0.00000 -4.39025 1 0.62832 0.41138 2 1.25664 -1.47331 3 1.88496 0.55312 4 2.51327 -0.69776 5 3.14159 0.00000 6 3.76991 0.00000 7 4.39823 1.47329 8 5.02655 -0.55312 9 5.65487 0.69776 10 6.28319 -5.00000 EXAMPLE N. 4 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 0.60000 1 0.60000 2 0.00000 3 -0.60000 4 -0.60000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 -1.00000 0.50000 0.60000 1 -0.86667 0.58000 0.60000 2 -0.73333 0.66000 0.60000 3 -0.60000 0.74000 0.60000 4 -0.46667 0.81997 0.59733 5 -0.33333 0.89630 0.53333 6 -0.20000 0.95840 0.38400 7 -0.06667 0.99490 0.14933 8 0.06667 0.99490 -0.14933 9 0.20000 0.95840 -0.38400 10 0.33333 0.89630 -0.53333 11 0.46667 0.81997 -0.59733 12 0.60000 0.74000 -0.60000 13 0.73333 0.66000 -0.60000 14 0.86667 0.58000 -0.60000 15 1.00000 0.50000 -0.60000 EXAMPLE N. 5 VALUES OF THE INPUT CONSTRAINTS I CONSTR(I) 0 1 1 1 2 2 3 1 4 1 VALUES OF THE OUTPUT CONSTRAINTS I CONSTR(I) 0 1 1 1 2 2 3 1 4 1 VALUES OF THE ERROR FLAGS : ERRC= 0 ; ERRE= 0 DIAGNOSTIC INFORMATION I DIAGN(I) 0 0 1 0 2 0 3 0 4 0 VALUES OF THE FIRST DERIVATIVES AT THE KNOTS I D(I) 0 1.00000 1 1.00000 2 0.50000 3 -0.50000 4 -1.00000 5 -1.00000 VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION POINTS I XTAB(I) Y0TAB(I) Y1TAB(I) Y2TAB(I) 0 1.27000 1.29175 1 -0.50000 -0.50000 2 5.70000 -0.70000 3 3.50000 1.55469 SHAR_EOF fi # end of overwriting check if test -f 'data1' then echo shar: will not over-write existing file "'data1'" else cat << \SHAR_EOF > 'data1' 10 .000000000E+00 .200000000E+00 .628318500E+00 .387785300E+00 .125663700E+01 .115105700E+01 .188495600E+01 .751056500E+00 .251327400E+01 .787785200E+00 .314159300E+01 -.200000100E+00 .376991100E+01 -.387785300E+00 .439823000E+01 -.115105700E+01 .502654800E+01 -.751056500E+00 .565486700E+01 -.787784900E+00 .628318500E+01 .200000200E+00 SHAR_EOF fi # end of overwriting check if test -f 'data2' then echo shar: will not over-write existing file "'data2'" else cat << \SHAR_EOF > 'data2' 8 0.0 0.0 1.0 0.75 2.0 1.0 3.0 0.75 4.0 0.0 5.0 -0.75 6.0 -1.0 7.0 -0.75 8.0 0.0 SHAR_EOF fi # end of overwriting check if test -f 'driver.f' then echo shar: will not over-write existing file "'driver.f'" else cat << \SHAR_EOF > 'driver.f' PROGRAM SBVSDR C ABSTRACT C C SBVSDR is a driver program which shows the use of the subroutines C SBVSIS, SBVSSC and SBVSSE with some examples. C C C DESCRIPTION C C SBVSDR consists of five examples which handle the execution of C SBVSSC (Boundary-Valued Shape-preserving Spline Computation) and C SBVSSE (Boundary-Valued Shape-preserving Spline Evaluation) or, when C convenient, the support routine SBVSIS, which simply collects the C calls to the previous ones, and cover the main features of the C package. The outputs are written to the external file RES. C For an easier reading of the program, all the output statements C are collected in the subroutine SWRTO. C The first and third examples make use of the data file data1 the C second of data2; obviously, they must be stored in the same C working directory as the executable version of SBVSDR. C The second and fourth examples make also use of external functions C which characterize the boundary conditions; although when the C non-separable boundary conditions are not imposed this function are C meaningless, vacuous function subprograms must be supplied for C linking purposes (for details see the description of the method and C of the input parameters BETA, BETAI, RHO and RHOI in SBVSIS). C C C SBVSDR makes use of: C C a) variables and arrays used as actual arguments for SBVSIS. C They have the same names of the corresponding dummy arguments C and their meaning is explained in the above subroutine; C C b) symbolic constant names, used for dimensioning the arrays. They C have the mnemonic form P-name, where P stands for PARAMETER and C 'name' is the name of a constant or variable which is used in C SBVSIS. So, their meaning can be easily obtained from the related C constant or variable 'name' described in the above subroutine; C C c) symbolic constants and variable names, used for input operations. C They will be explained at their first occurrence in executable C statements; C C d) function names used as actual arguments for SBVSIS. They have one C of the following forms SBn , SBnI , SRn , SRnI or SDUMFN. The C number n refers to the example in which the actual function is C used, B and R refer to functions 'beta' or 'rho' described in C SBVSIS and I , when present, to their inverse. SDUMFN is finally C used for linking purposes when 'beta' and/or 'rho' functions are C not referenced. C C The following are in class a) C C X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,KMAX,MAXSTP,XTAB,NTAB, C SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,Y0TAB,Y1TAB,Y2TAB,WORK,NWORK; C C in class b) C C PNP,PNTAB,PN,PCOMM,PPART,PNWORK; C C in class c) C C UN1,UN2,FNM1,FNM2; C C and in class d) C C SB2,SB2I,SR2,SR2I,SB4,SB4I,SDUMFN. C C SBVSDR calls the subroutines SBVSSC, SBVSSE and SWRTO C (at the entry points SWRT1A, SWRT1B, SWRT2, SWRT3A, SWRT3B, SWRT4, C SWRT5A, SWRT5B ). INTEGER PNP,PNTAB,PPART,PCOMM,PN,PNWORK PARAMETER (PNP=10,PNTAB=15,PPART=2,PCOMM=5,PN=6, * PNWORK=PCOMM+(PPART+7)*PNP+PN*(PN+11)/2+9) C The following statements define the symbolic names of constants used C for the input operations. INTEGER UN1,UN2 CHARACTER FNM1*9,FNM2*9 PARAMETER (UN1=1,UN2=2,FNM1='data1',FNM2='data2') INTEGER NP,N,K,OPT,CONSTR(0:PNP-1),KMAX,MAXSTP,NTAB,SBOPT, * Y0OPT,Y1OPT,Y2OPT,ERRE,ERRC,DIAGN(0:PNP-1),NWORK,I REAL X(0:PNP),Y(0:PNP),D0,DNP,D20,D2NP,EPS, * XTAB(0:PNTAB),D(0:PNP),D2(0:PNP), * Y0TAB(0:PNTAB),Y1TAB(0:PNTAB), * Y2TAB(0:PNTAB),WORK(1:PNWORK) REAL SB2,SB2I,SR2,SR2I,SB4,SB4I,SDUMFN EXTERNAL SB2,SB2I,SR2,SR2I,SB4,SB4I,SDUMFN C The following variables will be used in the output statements of C subroutine SWRTO. COMMON /SBVCOM/ NP,NTAB,ERRC,ERRE,X,Y,XTAB,Y0TAB,Y1TAB,Y2TAB, * D,D2,DIAGN,CONSTR C --------------------------------------------------------------------- C Call subroutine SWRTO for printing the heading in the output file. CALL SWRTO C --------------------------------------------------------------------- C Example n.1 . C C We want to construct a spline interpolating the set of data points C (x(i),y(i)) , i=0,1,...,10 , with x(i) = 2*i/20 and C y(i)= sin x(i) + ( cos 5x(i) )/5 , stored in the external file FNM1. C No constraints or boundary conditions are imposed on it and its C slopes are computed as the best approximation to the Bessel C estimates. The spline has continuity 1 and degree 3. Note that with C these selections we get the classical cubic Bessel interpolant. C We also want to estimate the spline at the set of tabulation points C xtab(i)=x(i) , i=0,...,10 . Then we want to evaluate the first and C the second derivative at the point 1.39 . C C C Read the data from the input file. It contains: C - the number NP of the data points, C - the coordinates (x(i),y(i)) of the i-th point, i=0,...,NP. OPEN (UNIT=UN1,FILE=FNM1,STATUS='OLD',ACCESS='SEQUENTIAL', * FORM='FORMATTED') READ (UNIT=UN1,FMT=*) NP READ (UNIT=UN1,FMT=*) (X(I),Y(I) , I=0,NP) CLOSE (UNIT=UN1,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=3 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: slopes as the best approximation to C the Bessel estimates, no boundary conditions, no constraints. OPT=210 C Set EPS, KMAX and MAXSTP to 0 for the automatic choice of their C values. EPS=0.E+00 KMAX=0 MAXSTP=0 C Call SBVSSC to compute the spline's parameters. CALL SBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,SDUMFN, * SDUMFN,SDUMFN,SDUMFN,KMAX,MAXSTP,ERRC,D,D2,DIAGN, * WORK,NWORK) C Set the options for evaluating the spline at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=10 DO 10 I=0,NTAB XTAB(I)=X(I) 10 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call SBVSSE for evaluating the spline. CALL SBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call SWRTO at the entry point SWRT1A for printing the results in the C output file. CALL SWRT1A C Set the options for evaluating the spline at the tabulation point. Y0OPT=0 Y1OPT=1 Y2OPT=1 C Define the new tabulation point. NTAB=0 XTAB(0)=0.139E+01 C Set SBOPT to 2 for selecting the binary search algorithm. SBOPT=2 C Call SBVSSE for evaluating the spline. CALL SBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call SWRTO at the entry point SWRT1B for printing the results in the C output file. CALL SWRT1B C --------------------------------------------------------------------- C Example n.2 . C C Given the set of points (x(i),y(i)) , i=0,...,8 where C {x(0),...,x(8)}={0.,1.,2.,3.,4.,5.,6.,7.,8.} and C {y(0),...,y(8)}={0.,0.75,1.,0.75,0.,-0.75,-1.,-0.75,0.} C stored in the external file FNM2, we want to construct a C co-monotone spline subject to the additional boundary condition C s'(8)=s'(0) , s''(8)=s''(0) . C We also ask the slopes to be the best approximation to the set of C third order monotonicity-preserving estimates described in SBVSSC C (see the input parameter OPT and the subroutine TODC). C The spline has degree 6 and continuity 2. C We evaluate the spline and its second derivative at the set of C tabulation points xtab(i), i=0,...,10 uniformly distributed in the C interval (x(0),x(8)). C We use in this case the support routine SBVSISP. C C Read the data from the input file. It contains: C - the number NP of the data points, C - the coordinates (x(i),y(i)) of the i-th point, i=0,...,8 . OPEN (UNIT=UN2,FILE=FNM2,STATUS='OLD', * ACCESS='SEQUENTIAL',FORM='FORMATTED') READ (UNIT=UN2,FMT=*) NP READ (UNIT=UN2,FMT=*) (X(I),Y(I),I=0,NP) CLOSE (UNIT=UN2,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=6 K=2 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: slopes as the best approximation to C the set of third order accurate estimates, non-separable boundary C conditions, monotonicity constraints. OPT=321 C Set EPS, KMAX and MAXSTP to 0 for the automatic choice of their C values. EPS=0.0E+00 KMAX=0 MAXSTP=0 C Set the options for evaluating the spline and its second derivative C at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=1 C Define the number and the set of tabulation points. NTAB=10 DO 20 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 20 CONTINUE C Set SBOPT to 1 for the sequential search algorithm. SBOPT=1 C Call SBVSIS for computing the spline. CALL SBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS, * SB2,SB2I,SR2,SR2I,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call SWRTO at the entry point SWRT2 for printing the results into C the output file. CALL SWRT2 C --------------------------------------------------------------------- C Example n.3 . C C Given the set of data points (x(i),y(i)) i=0,1,...,10 as in the C example 1, we want to construct a spline of degree 6 and continuity 2 C which is co-convex and subject to separable boundary conditions: C s'(x(0))=y'(x(0)) , s'(x(10))=y'(x(10)) , s''(x(0))=y''(x(0)) , C s''(x(10))=y''(x(10)) . The parameters are determined only by the C convexity constraints. We evaluate the spline and its first C derivative at the tabulation points xtab(i) , i=0,...,15 uniformly C distributed in the interval (x(0),x(10)) and then the second C derivative at the knots. C C C Read the data from the input file. OPEN (UNIT=UN1,FILE=FNM1,STATUS='OLD', * ACCESS='SEQUENTIAL',FORM='FORMATTED') READ (UNIT=UN1,FMT=*) NP READ (UNIT=UN1,FMT=*) (X(I),Y(I),I=0,NP) CLOSE (UNIT=UN1,STATUS='KEEP') C Set the degree and the class of continuity of the spline. N=6 K=2 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: the slopes are computed only from C the constraints, non-separable boundary conditions, C convexity constraints. OPT=132 C Set the boundary conditions. D0=COS(X(0))+SIN(5.0*X(0)) DNP=COS(X(NP))+SIN(5.0*X(NP)) D20=-SIN(X(0))-5.0*COS(5.0*X(0)) D2NP=-SIN(X(NP))-5.0*COS(5.0*X(NP)) C Set EPS, KMAX and MAXSTP for the automatic choice of their values. EPS=0.0E+00 KMAX=0 MAXSTP=0 C Call SBVSSC for computing the spline. CALL SBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,SDUMFN, * SDUMFN,SDUMFN,SDUMFN,KMAX,MAXSTP,ERRC,D,D2,DIAGN, * WORK,NWORK) C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=1 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=15 DO 30 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 30 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call SBVSSE for evaluating the spline. CALL SBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call SWRTO at the entry point SWRT3A for printing the results in the C output file. CALL SWRT3A C Set the options for evaluating the second derivative at the knots. Y0OPT=0 Y1OPT=0 Y2OPT=1 C Define the number and the set of tabulation points. NTAB=NP DO 40 I=0,NTAB XTAB(I)=X(I) 40 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call SBVSSE for evaluating the spline. CALL SBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C Call SWRTO at the entry point SWRT3B for printing the results in the C output file. CALL SWRT3B C --------------------------------------------------------------------- C Example n.4 . C C Given the 'Runge function': f(x) = 1/(1+x**2) , we want to construct C a co-monotone and co-convex spline which interpolates the function at C x(i) = -1+2*i/4 , i=0,1,...,4 . We also impose the additional C boundary conditions s'(1)=-s'(-1) and we compute the slopes as the C best approximation to the estimates s'(x(i)) = f'(x(i)) , i=0,...,4. C The spline has degree 3 and continuity 1. We evaluate the spline and C its first derivative at the tabulation points xtab(i), i=0,...,15 C uniformly distributed in the interval(x(0),x(4)). C It is used the support routine SBVSIS C C Compute the input data NP=4 DO 50 I=0,NP X(I)=-1.0+2.0*I/4.0 Y(I)=1.0/(1.0+X(I)**2) D(I)=-2.0*X(I)/((1.0+X(I)**2)**2) 50 CONTINUE C Set the degree and the class of continuity of the spline. N=3 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: best approximation to user supplied C derivatives, boundary conditions, monotonicity and convexity C constraints. OPT=423 C Set EPS, KMAX and MAXSTP EPS=1.0E-05 KMAX=5 MAXSTP=15 C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=1 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=15 DO 60 I=0,NTAB XTAB(I)=X(0)+I*(X(NP)-X(0))/NTAB 60 CONTINUE C Set SBOPT to 1 for selecting the sequential search algorithm. SBOPT=1 C Call SBVSIS for computing the spline. CALL SBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,SB4,SB4I, * SDUMFN,SDUMFN,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call SWRTO at the entry point SWRT4 for printing the results into C the output file. CALL SWRT4 C --------------------------------------------------------------------- C Example n.5 . C C Given the points { (0,0),(1,1),(2,2),(3,2),(4,1),(5,0) } we want to C construct a spline of degree 4 and continuity 1 which is co-convex C in the interval (2,3) and co-monotone elsewhere. We also require the C separable boundary conditions s'(0)=1 , s'(5)=-1. The derivatives C are computed as the best approximation to Bessel estimates. C The spline is evaluated at four random points. C The routine SBVSIS is used. C C Compute the input data NP=5 X(0)=0.0E+00 X(1)=1.0E+00 X(2)=2.0E+00 X(3)=3.0E+00 X(4)=4.0E+00 X(5)=5.0E+00 Y(0)=0.0E+00 Y(1)=1.0E+00 Y(2)=2.0E+00 Y(3)=2.0E+00 Y(4)=1.0E+00 Y(5)=0.0E+00 C Set the degree and the class of continuity of the spline. N=4 K=1 C Set the size of the work area NWORK=PCOMM+(PPART+7)*NP+(N*(N+11))/2+9 C Set the options on the spline: Bessel like derivatives, C separable boundary conditions, user specified shape constraints. OPT=234 C Set the boundary conditions. D0=1.0E+00 DNP=-1.0E+00 C Select the shape constraints. CONSTR(0)=1 CONSTR(1)=1 CONSTR(2)=2 CONSTR(3)=1 CONSTR(4)=1 C Call SWRTO at the entry point SWRT5A for printing the input C constraints values. CALL SWRT5A C Set EPS, KMAX and MAXSTP for the automatic choice of their values. EPS=0. KMAX=0 MAXSTP=0 C Set the options for evaluating the spline and its first derivative C at the tabulation points. Y0OPT=1 Y1OPT=0 Y2OPT=0 C Define the number and the set of tabulation points. NTAB=3 XTAB(0)=0.127E+01 XTAB(1)=-0.5E+00 XTAB(2)=0.57E+01 XTAB(3)=0.35E+01 C Set SBOPT to 2 for selecting the binary search algorithm. SBOPT=2 C Call SBVSIS for computing the spline. CALL SBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,SDUMFN, * SDUMFN,SDUMFN,SDUMFN,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2,DIAGN, * Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C Call SWRTO at the entry point SWRT5B for printing the results into C the output file. CALL SWRT5B STOP END C --------------------------------------------------------------------- SUBROUTINE SWRTO C SWRTO is an auxiliary routine which prints the results of the driver C program SBVSDR in the external file RES. All the parameters C are shared with the main program in the common area and their meaning C is explained in SBVSDR. C C Items of possible interest are: C C OUN : symbolic name of an integer constant, used for identifying C the logical output number; C FNAME : symbolic name of a character constant, used for identifying C the name of the output file. C C The subroutine is called at the entry points SWRTO, SWRT1A, SWRT1B, C SWRT2, SWRT3A, SWRT3B, SWRT4, SWRT5A, SWRT5B. INTEGER PNP,PNTAB,PN PARAMETER (PNP=10,PNTAB=15,PN=6) INTEGER OUN CHARACTER FNAME*12 PARAMETER (OUN=10,FNAME='RES') INTEGER NP,NTAB,CONSTR(0:PNP-1),ERRC,DIAGN(0:PNP-1),ERRE INTEGER I REAL X(0:PNP),Y(0:PNP),XTAB(0:PNTAB),Y0TAB(0:PNTAB), * Y1TAB(0:PNTAB),Y2TAB(0:PNTAB),D(0:PNP), * D2(0:PNP) COMMON /SBVCOM/ NP,NTAB,ERRC,ERRE,X,Y,XTAB,Y0TAB,Y1TAB,Y2TAB, * D,D2,DIAGN,CONSTR C Connect the output file FNAME to the logical unit number OUN. OPEN (UNIT=OUN,FILE=FNAME,STATUS='UNKNOWN', * ACCESS='SEQUENTIAL',FORM='FORMATTED') C Write the heading of the output file WRITE(UNIT=OUN,FMT=10) RETURN C Write the results of the example 1. ENTRY SWRT1A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 1' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=100) (I,XTAB(I),Y0TAB(I),I=0,NTAB) RETURN ENTRY SWRT1B WRITE (UNIT=OUN,FMT=25) 'FURTHER EVALUATIONS' WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=011) (I,XTAB(I),Y1TAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 2. ENTRY SWRT2 WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 2' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=45) WRITE (UNIT=OUN,FMT=55) (I,D(I),D2(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=101) (I,XTAB(I),Y0TAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 3. ENTRY SWRT3A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 3' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=45) WRITE (UNIT=OUN,FMT=55) (I,D(I),D2(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=110) (I,XTAB(I),Y0TAB(I),Y1TAB(I),I=0,NTAB) RETURN ENTRY SWRT3B WRITE (UNIT=OUN,FMT=25) 'FURTHER EVALUATIONS' WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=001) (I,XTAB(I),Y2TAB(I),I=0,NTAB) RETURN C Write the results of the example 4. ENTRY SWRT4 WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 4' WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=110) (I,XTAB(I),Y0TAB(I),Y1TAB(I),I=0,NTAB) RETURN C Write the results of the example 5. ENTRY SWRT5A WRITE (UNIT=OUN,FMT=20) 'EXAMPLE N. 5' WRITE (UNIT=OUN,FMT=70) WRITE (UNIT=OUN,FMT=90) (I,CONSTR(I),I=0,NP-1) RETURN ENTRY SWRT5B WRITE (UNIT=OUN,FMT=80) WRITE (UNIT=OUN,FMT=90) (I,CONSTR(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=30) ERRC,ERRE WRITE (UNIT=OUN,FMT=33) WRITE (UNIT=OUN,FMT=36) (I,DIAGN(I),I=0,NP-1) WRITE (UNIT=OUN,FMT=40) WRITE (UNIT=OUN,FMT=50) (I,D(I),I=0,NP) WRITE (UNIT=OUN,FMT=60) WRITE (UNIT=OUN,FMT=100) (I,XTAB(I),Y0TAB(I),I=0,NTAB) C Disconnect the output file. CLOSE (UNIT=OUN,STATUS='KEEP') RETURN 10 FORMAT(//22X,'SRESULTS.OUT'// * ' THIS FILE CONTAINS THE COMPUTED RESULTS OF THE DRIVER'/ * ' PROGRAM SBVSDR, WHICH DEMONSTRATES THE USE OF THE'/ * ' SUBROUTINES SBVSIS, SBVSSC AND SBVSSE WITH SOME EXAMPLES') 20 FORMAT(//22X,A13/) 25 FORMAT(/22X,A19/) 30 FORMAT(/' VALUES OF THE ERROR FLAGS : ERRC=', * I2,' ; ERRE=',I2) 33 FORMAT(/' DIAGNOSTIC INFORMATION'//, * 10X,'I',7X,'DIAGN(I)'/) 36 FORMAT(9X,I2,11X,I1) 40 FORMAT(/' VALUES OF THE FIRST DERIVATIVES AT THE KNOTS'//, * 10X,'I',9X,'D(I)'/) 45 FORMAT(/' VALUES OF THE FIRST AND SECOND DERIVATIVES AT ', * 'THE KNOTS'//,10X,'I',9X,'D(I)',9X,'D2(I)'/) 50 FORMAT(9X,I2,5X,F9.5) 55 FORMAT(9X,I2,5X,F9.5,5X,F9.5) 60 FORMAT(/ * ' VALUES OF THE SPLINE AND ITS DERIVATIVES AT THE TABULATION', * ' POINTS'//, * 7X,'I',5X,'XTAB(I)',4X,'Y0TAB(I)',4X,'Y1TAB(I)',4X,'Y2TAB(I)'/) 70 FORMAT(/' VALUES OF THE INPUT CONSTRAINTS'//, * 10X,'I',7X,'CONSTR(I)'/) 80 FORMAT(/' VALUES OF THE OUTPUT CONSTRAINTS'//, * 10X,'I',7X,'CONSTR(I)'/) 90 FORMAT(9X,I2,11X,I1) 001 FORMAT(5X,I3,3X,F9.5,27X,F9.5) 011 FORMAT(5X,I3,3X,F9.5,15X,F9.5,3X,F9.5) 100 FORMAT(5X,I3,3X,F9.5,3X,F9.5) 101 FORMAT(5X,I3,3X,F9.5,3X,F9.5,15X,F9.5) 110 FORMAT(5X,I3,3X,F9.5,3X,F9.5,3X,F9.5) END C --------------------------------------------------------------------- C The following functions define the non-separable boundary conditions. C For linking purposes, these must be given even when the conditions C are not required. C In this case the 'vacuous function' SDUMFN is used. REAL FUNCTION SB2(X) REAL X SB2=X RETURN END REAL FUNCTION SB2I(Y) REAL Y SB2I=Y RETURN END REAL FUNCTION SR2(X) REAL X SR2=X RETURN END REAL FUNCTION SR2I(Y) REAL Y SR2I=Y RETURN END REAL FUNCTION SB4(X) REAL X SB4=-X RETURN END REAL FUNCTION SB4I(Y) REAL Y SB4I=-Y RETURN END REAL FUNCTION SDUMFN(T) REAL T SDUMFN=0.0E+00*T RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' SUBROUTINE DBVSIS (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS, * BETA,BETAI,RHO,RHOI,KMAX,MAXSTP,XTAB,NTAB, * SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,D,D2, * DIAGN,Y0TAB,Y1TAB,Y2TAB,WORK,NWORK) C DBVSIS is merely a support routine which calls DBVSSC and DBVSSE C for the computation of the needed parameters and for the evaluation C of a shape-preserving, C(k), k=1,2 , interpolating spline, C optionally subject to boundary conditions. C The use of DBVSIS is not recommended when more evaluations of the C same spline are required; in this case it is better to separately C call DBVSSC and then DBVSSE repeatedly. C For an explanation of input and output parameters, the user is C referred to the comments of DBVSSC and DBVSSE. EXTERNAL BETA,BETAI,RHO,RHOI INTEGER NP,N,K,OPT,CONSTR(0:NP-1),NTAB,KMAX,MAXSTP,SBOPT, * Y0OPT,Y1OPT,Y2OPT,DIAGN(0:NP-1),ERRC,ERRE,NWORK DOUBLE PRECISION X(0:NP),Y(0:NP),D0,DNP,D20,D2NP,XTAB(0:NTAB), * EPS,BETA,BETAI,RHO,RHOI,D(0:NP),D2(0:NP), * Y0TAB(0:NTAB),Y1TAB(0:NTAB),Y2TAB(0:NTAB), * WORK(1:NWORK) CALL DBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,BETA, * BETAI,RHO,RHOI,KMAX,MAXSTP,ERRC,D,D2,DIAGN, * WORK,NWORK) CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) RETURN END C --------------------------------------------------------------------- SUBROUTINE DBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS, * BETA,BETAI,RHO,RHOI,KMAX,MAXSTP,ERRC,D,D2, * DIAGN,WORK,NWORK) C ------------------------------------------------- C Lines 49-549 are comment lines. C The actual code begins at line 555. C ------------------------------------------------- C ABSTRACT: C C DBVSSC is designed to compute the coefficients (first and, if C appropriate, second derivatives) of a shape-preserving spline, of C continuity class C(k), k=1,2 , which interpolates a set of data C points and, if required, satisfies additional boundary conditions. C DBVSSC furnishes the input parameters for DBVSSE, which, in turn, C provides to evaluate the spline and its derivatives at a set of C tabulation points. C C The user is allowed to use the following options: C C - to compute a spline subject to: C - no constraint, C - monotonicity constraints, C - convexity constraints, C - monotonicity and convexity constraints, C - one of the above constraints in each subinterval; C C - to impose separable or non-separable boundary conditions on the C spline, C C - to assign the first derivatives d(i), i=0,1,...,np , in input or to C compute them from the constraints only or as the best approximation C to a set of optimal values. Although the final sequence of C derivatives does in any case satisfy the imposed restrictions on C the shape, the resulting graphs may exhibit different behaviors. C C C REMARK: C C In these comments variable and array names will be denoted with C capital letters, and their contents with small letters. Moreover: C .IN. means belonging to; C .INT. stands for intersection. C C C The code has the following structure: C C DBVSSC C DBVC C DSTINF C DMSK1 C DMSK2 C DPRJ0 C DPRJ1 C DTDC C DMNMOD C DMDIAN C DALG3 C DPRJ0 C DALG1 C DPRJ1 C DINTRS C DTST C DFPSVF C DSL C DALG1D C DPRJ1 C DAL2 C DPRJ2 C DAL2DP C DMNIND C DPRJ2 C DSL C DSCDRC C C C CALLING SEQUENCE: C C CALL DBVSSC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,EPS,BETA, C * BETAI,RHO,RHOI,KMAX,MAXSTP,ERRC,D,D2,DIAGN, C * WORK,NWORK) C C C INPUT PARAMETERS: C C X : floating array, of bounds 0:NP, containing the data C abscissas x(i), i=0,1,...,np. C Restriction: x(i).LT.x(i+1), i=0,1,...,np. C Y : floating array, of bounds 0:NP, containing the data C ordinates y(i), i=0,1,...,np. C NP : integer variable, defining the number of interpolation C points. Restriction: np.GE.2 . C N : integer variable, containing the degree of s. C Restriction: n.GE.3 . C K : integer variable, containing the class of continuity of s. C Restriction: k.EQ.1 or k.EQ.2 and n.GE.3*k . C OPT : integer variable, containing a control parameter. It is C a three-digit decimal of the form pqr (that is of C numerical value p*100+q*10+r ) where: C r controls the constraints on the shape. C q controls the boundary conditions and C p controls the computation of the derivatives, C More specifically: C r=0 (opt=pq0) : no constraint on the shape is required; C r=1 (opt=pq1) : monotonicity constraints are required; C r=2 (opt=pq2) : convexity constraints are required; C r=3 (opt=pq3) : monotonicity and convexity constraints are C required; C r=4 (opt=pq4) : local constraints for any subinterval are C supplied by the user (see the description C of the array CONSTR); C q=1 (opt=p1r) : no boundary condition is imposed; C q=2 (opt=p2r) : non-separable boundary conditions are C imposed (see the description of BETA, C BETAI, RHO, RHOI); C q=3 (opt=p3r) : separable boundary conditions are imposed C (see the description of D0, DNP, D20, C D2NP); C p=1 (opt=1qr) : the sequence of first derivatives C d(0),....,d(np) is computed using the C constraints only using subroutine DAL2; C p=2 (opt=2qr) : the sequence is computed as the constrained C best approximation to Bessel derivatives C using subroutine DAL2DP; C p=3 (opt=3qr) : the sequence is computed as the constrained C best approximation to a set of third order C accurate derivative estimates produced in C subroutine DTDC using subroutine DAL2DP C (since this estimates are inherently mono- C tonicity preserving, it is not recommended C to associate this option with the convexity C constraints only); C p=4 (opt=4qr) : the sequence is computed as the constrained C best approximation to a set of values given C in input by the user using DAL2DP; note C that opt.EQ.410 will provide the classical C C(k) function interpolating the data and C the derivatives. C Restriction: ( p.GE.1 .AND. p.LE.4 ) .AND. C ( q.GE.1.AND. q.LE.3 ) .AND. C ( r.GE.0 .AND. r.LE.4 ) .AND. C .NOT. ( r.EQ.0 .AND. p.EQ.1 ) . C D0 : floating variable containing the left separable boundary C condition for the first derivative (d(0)=d0). C D0 is referenced only when q=3 . C DNP : floating variable containing the right separable boundary C condition for the first derivative (d(np)=dnp). C DNP is referenced only when q=3 . C D20 : floating variable containing the left separable boundary C condition for the second derivative (d2(0)=d20). C D20 is referenced only when q=3 and k=2 . C D2NP : floating variable containing the right separable boundary C condition for the second derivative (d2(np)=d2np). C D2NP is referenced only when q=3 and k=2 . C EPS : floating variable, containing a value for computing the C relative tolerance of the method. It should be set greater C or equal to the machine precision. However, if eps.LE.0, C DBVSSC resets it to 0.0001 which has turned out to be a C good choice for graphical applications. C BETA : user supplied function, which represents non-separable C boundary conditions for the first derivatives. C BETA is referenced only when q=2 . C BETAI : user supplied function, which is the inverse of BETA. C BETAI is referenced only when q=2 . C RHO : user supplied function, which represents non-separable C boundary conditions for the second derivatives. C RHO is referenced only when q=2 and k=2 . C RHOI : user supplied function, which is the inverse of RHO. C RHOI is referenced only when q=2 and k=2 . C KMAX : integer variable, containing the number of iterations C allowed for selecting the minimal set ASTAR described C below. If kmax.LE.0, DBVSSC resets it to 10 . C KMAX is referenced only when q=2 . C MAXSTP : integer variable, containing the number of steps allowed C to find dstar in the set of admissible values. C If maxstp.LE.0, DBVSSC resets it to 10 . C MAXSTP is referenced only when q=2 . C C C INPUT / OUTPUT PARAMETERS: C C CONSTR : integer array, of bounds 0:NP , containing, in input the C desired constraints on the shape for each subinterval. C constr(i)=kind , kind=0,1,2,3 , means that none, monotoni- C city, convexity, monotonicity and convexity constraint is C imposed on the subinterval [x(i),x(i+1)]. If constr(i) is C not compatible with the data it is relaxed according to C their shape (see subroutine DMSK1 for details). So, on out- C put, CONSTR contains the constraints actually imposed. C For example, if the data are convex and on input we have C constr(i)=3 , the result in output will be constr(i)=2. C Restriction: constr(i).GE.0 .AND. constr(i).LE.3 . C CONSTR is referenced only when r=4 . C D : floating array, of bounds 0:NP, containing the first C derivatives at x(i), i=0,1,...,np . If p=4 , d(i) is the C input value to be approximated by the computed derivative, C which is then stored in the same location. C On output, D is computed by the routine DAL2 if p=1 and C is computed by the routine DAL2DP if p=2 or p=3 . C C C OUTPUT PARAMETERS C C ERRC : integer variable, containing an error flag which displays C the status of the code. The status is divided into: severe C error (error on the input data, no computation has been C done), error (some computation has been done and some C information or suggestions are available), warning (some C requirement is not fulfilled, but the spline's parameters C have been computed), success. C errc=0 : success, normal return of the code; C errc=1 : severe error, incorrect assignment for some of C the values nwork, opt, np; C errc=2 : severe error, for some i the restriction C 0.LE.constr(i) .AND. constr(i).LE.3 is not C fulfilled; C errc=3 : severe error, incorrect assignment for some of C the values n,k; C errc=4 : severe error, the restriction x(i).LT.x(i+1) is C not fulfilled for some i; C errc=5 : error, the problem does not have any solution C because the set C betai ( phi(a(0,k)) .INT. beta(a(0,k)) ) C is empty for some k. In other words the boundary C conditions cannot be satisfied and the output C parameters are meaningless. C The user is suggested to increase the value of n. C errc=6 : warning; for some i, the constraints on the C interval [x(i),x(i+1)] are too strong and they C have not been considered. There is no guarantee C that the spline is shape-preserving within all C the intervals. More accurate diagnostic details C can be found in the array DIAGN. C The user is suggested to increase the value of n. C errc=7 : error, dstar such that beta(dstar).IN.phi(dstar) C has not been found. The integer parameter maxstp C should be increased. C The output parameters are meaningless. C errc=8 : error, both situations described in errc=6 and C errc=7 have occurred. C errc=9 : warning, one of the separable boundary conditions C d(0)=d0 and/or d(np)=dnp are not compatible C with the constraints in [x(0),x(1)] and/or C [x(np-1),x(np)] which have consequently been C relaxed. The user is suggested to increase the C value of n. More accurate diagnostic details can C be found in the array DIAGN. C errc=10: warning, both situations described for errc=6 and C errc=9 have occurred. C errc=11: warning, one of the separable boundary conditions C d2(0)=d20 and/or d2(np)=d2np are not compatible C with the constraints in [x(0),x(1)] and/or C [x(np-1),x(np)] . The boundary conditions have C consequently been approximated. The user is C suggested to increase the value of n. C errc=12: warning, both situations described for errc=6 and C errc=11 have occurred. C errc=13: warning, both situations described for errc=9 and C errc=11 have occurred. C errc=14: warning, both situations described for errc=10 and C errc=11 have occurred. C errc=15: warning, the non-separable boundary conditions C d2(np)=rho(d2(0)) are not compatible with the C constraints. The boundary conditions have C consequently been approximated. The user is C suggested to increase the value of n. C errc=16: warning, both situations described for errc=6 and C errc=15 have occurred. C errc=17: warning, both situations described for errc=9 and C errc=15 have occurred. C errc=18: warning, both situations described for errc=10 and C errc=15 have occurred. C D2 : floating array of bounds 0:NP containing the second C derivatives at knots. D2 is computed in subroutine DCDERC . C D2 is referenced only when k=2 . C DIAGN : integer array of bounds 0:NP-1 containing further C diagnostic information: C diagn(i)=0 : the constraints in the interval [x(i),x(i+1)] C have been satisfied; C diagn(i)=1 : the constraints in the interval [x(i),x(i+1)] C have not been satisfied; C C C C OTHER PARAMETERS: C C WORK : floating array, of bounds 1:NKORK, which is used as C a work area to store intermediate results. C The same array can be used to provide workspace for both C the main subroutines DBVSSC and DBVSSE . C NWORK : integer variable containing the size of the work area. C Restriction: nwork .GE. comm+(part+7)*np+(n*(n+11))/2+9 C that is C nwork .GE. 5+(2+7)*np+(n*(n+11))/2+9 C C C ------------------------------------------------ C C METHOD: C C Let the integers n and k, such that k=1,2 ; n >= 3k , and the C sequence of points (x(i),y(i)), i=0,1,...,np , with C x(0) < x(1) < ... = 0, v >= 0, v =< -u+ n/k delta(i) } C (3) C ( D(i) := { (u,v).IN.RxR : u =< 0, v =< 0, v >= -u+ n/k delta(i) } ) C C s is convex (CVX) ( concave (CNC) ) if and only if (d(i),d(i+1)) C belongs to C C D(i) := { (u,v).IN.RxR : v >= - (k/(n-k)) u + (n/(n-k)) delta(i) , C v =< - ((n-k)/k) u + (n/k) delta(i) } C (4) C ( D(i) := { (u,v).IN.RxR : v =< - (k/(n-k)) u + (n/(n-k)) delta(i) , C v >= - ((n-k)/k) u + (n/k) delta(i) } ) C C and that s is I (D) and CVX (CNC) if and only if (d(i),d(i+1)) C belongs to C C D(i) := { (u,v) satisfying (3) and (4) } . C C So, if we choose the family of sets D(i) , i=0,1,...,np-1 , according C to the shape of the data, we have to solve: C C PROBLEM P1. Does a sequence ( d(0), d(1), ..., d(np) ) such that C (d(i),d(i+1)) .IN. D(i) , i=0,1,...,np-1 , exist ? C C PROBLEM P2. If P1 is feasible, how can a (the best) solution be C computed ? C C Let DPRJ1: RxR -> R and DPRJ2: RxR -> R be, respectively, the C projection maps from uv-plane onto the u-axis and v-axis and let us C denote with B(i) := DPRJ1(D(i)) : C C ALGORITHM A1[B0]. C 1. Set A(0):=B(0); J:=np. C 2. For i=1,2,...,np C 2.1. Set A(i):= DPRJ2( D(i-1).INT.{ A(i-1) x B(i) } ) . C 2.2. If A(i) is empty, set J:=i and stop. C 3. Stop. C C We have the following result: C C THEOREM 1. P1 has a solution if, and only if, J=np, that is A(i) is C not empty , i=0,1,...,np . If ( d(0), d(1), ...,d(np) ) C is a solution then d(i).IN.A(i) , i=0,1,...,np . C C A solution can be computed with the following algorithm: C C ALGORITHM A2[A(np),B0]. C 1. Choose d(np).IN.A(np). C 2. For i=np-1, np-2, ..., 0 C 2.1. Choose d(i).IN.DPRJ1( D(i).INT.{ A(i) x { d(i+1) }}). C 3. Stop. C C For more theoretical details about A1 and A2 see \1\ , and for C practical details see subprograms DALG1, DAL2, DAL2DP. In the latter C a dynamic programming scheme is used to find the best solution in C the feasible set. More specifically, it is possible to compute the C values d(i),i=0,..,np which satisfy the constraints and are as close C as possible to another sequence which does not satisfy the C constraints but is, in some sense, optimal. C C From a theoretical point of view, algs A1 and A2 give a complete C answer to problems P1 and P2. However, it could be pointed out that, C for practical applications, we would like to have the best possible C plot, whether or not P1 has solutions. Let us suppose that the C problem is solvable from 0 to j and from j to np, but that alg A1 C applied to the whole family of sets D(i), i=0,1,...,np-1 gives C J.eq.j.ne.np ; if we reset D(j-1) := A(j-1) x B(j) , alg A1 applied C to this new family of sets will produce J=np . However, it must be C recalled that, in this way, we do not consider the constraints in the C interval [x(j-i),x(j)] and so there is no guarantee that the spline C is shape-preserving in this interval. Whenever this fact cannot be C accepted it is convenient to rerun the code with a larger value for C the degree n , because the domains of constraints enlarge as n C increases (see (3) and (4)). C C It is immediate to see that separable boundary conditions of the form C C d(0) := d0 ; d(np) := dnp C C can be easily inserted with a reduction of the corresponding C admissible sets which does not modify the above theory: C C D(0) := D(0).INT.{d(0)=d0} ; D(np) := D(np).INT.{d(np)=dnp} C C In the case k=2 the corresponding conditions d2(0) = d20 , C d2(np) = d2np are imposed only if not in contrast with the shape of C the data; otherwise the admissible values for d2(0) and d2(np) C respectively closest to d20 and d2np are chosen. C C Now, let beta be a continuous function from R to R, with continuous C inverse betai, we want to solve the following non-separable boundary C valued problem: C C PROBLEM P3. Do sequences ( d(0), d(1), ..., d(np) ) , such that C (d(i),d(i+1)).IN.D(i), i=0,1,...,np-1 and C d(np) = beta ( d(0) ) , exist ? C C It is obvious that a solution of this new problem, if it exists, can C be found amongst the solutions of P1. Let A(0), A(1),...,A(np) be the C sequence of sets given by alg A1 (we assume that A(i) is not empty, C i=0,1,...,np , that is P1 is solvable or, if this is not the case, C the constraints have been relaxed ), we can assume that C A(np) = phi(A(0)) , where phi: R -> R is a set valued function C (see \1\ for details). It can be demonstrated that: C C THEOREM 2. P1 is solvable if, and only if, there is dstar.IN.A(0) C such that beta(dstar).IN.phi({dstar}) . C C It should be noted that if ( d(0), d(1), ..., d(np) ) satisfies P1, C d(0) .IN. betai(phi(A(0)).INT.beta(A(0))) =: gamma(A(0)) C and, consequently, the set of admissible values is reduced. If we C repeat this procedure, we get a gradually diminishing admissible set C for d(0). We define C ASTAR := lim A(0,m) where C A(0,0) := A(0) and A(0,m) := gamma(A(0,m-1)) ; C ASTAR is the minimal admissible set for dstar. We can now combine the C various theorems and algorithms and give the general algorithm to C solve P3: C C ALGORITHM A3. C 1. Set A(0,0) := B0 ; m:=1. C 2. Use A1[A(0,0)] for computing phi (A(0,0)). C 3. Set A(0,1) := gamma(A(0,0)) C = betai(phi(A(0,0)).INT.beta(A(0,0))). C 4. If A(0,1) is empty, stop (P1 is unsolvable). C 5. While ( convergence test not satisfied ) do C 5.1. Use A1[A(0,m)] for computing A(np,m) = phi (A(0,m)). C 5.2. Set A(0,m+1) := gamma(A(0,m)). C 5.3. Set m:=m+1. C 6. Set ASTAR := A(0,m). C 7. Use A1[{d(0)}] to find dstar.IN.ASTAR such that C beta(dstar).IN.phi(dstar). C 8. Use A2[beta(dstar),dstar] for computing a sequence C ( d(0), d(1), ..., d(np) ) which solves P1. C C In the case k=2 the corresponding condition d2(np) = beta2(d2(0)) C is imposed only if not in contrast with the shape of C the data; otherwise the admissible values for d2(0) and d2(np) C closest to the boundary condition are chosen. C C References C C \1\ P.Costantini: Boundary Valued Shape-Preserving Interpolating C Splines, ACM Trans. on Math. Softw., companion paper. C \2\ R.Bellman, S.Dreyfus: Applied Dynamic Programming, Princeton C University Press, New York, 1962. C \3\ H.T.Huynh: Accurate Monotone Cubic Interpolation, SIAM J. Num. C Anal., 30 (1993), 57-100. C C The ideas involved in Algorithm A3 have been implemented in the code C in a general form. Since Algorithm A3 resembles closely the abstract C formulation it could, therefore, be used for several practical C problems. The particular case actually treated is reflected in the C contents of the information array INFO (see its description in C subroutine DSTINF) which contains all the data needed for the C construction of the operators DPRJ0, DPRJ1 and DPRJ2. C C As a consequence, the user has the following options: C C - to compute a Spline subject to: C - no constraint; C - monotonicity constraints, C - convexity constraints, C - monotonicity and convexity constraints, C - one of the above constraints in each subinterval, as C specified in the corresponding array CONSTR; C C - to impose separable or non-separable boundary conditions on the C spline. In the latter case, the external functions BETA, BETAI, C RHO and RHOI must be supplied, C C - to assign the first derivatives d(i), i=0,1,...,np , in input or to C compute them from the only constraints or as the best approximation C to a set of optimal values. Although the final sequence of C derivatives does in any case satisfy the imposed restrictions on C the shape, the resulting graphs may exhibit different behaviors. C C See the description of the input parameter OPT for more details. C ------------------------------------------------ C End of comments. C ------------------------------------------------ INTEGER COMM,PART C COMM contains the number of global data referring to the initial C points (x(i),y(i)) stored in the array INFO, described in C subroutine DSTINF. C PART contains the number of particular data referring to each C interval (x(i),x(i+1)) , i=0,1,...,np , stored in the array INFO. PARAMETER (COMM=5, PART=2) EXTERNAL BETA,BETAI,RHO,RHOI INTEGER NP,N,K,OPT,CONSTR(0:NP-1),KMAX,MAXSTP,ERRC, * DIAGN(0:NP-1),NWORK,I1,I2,I3,I4,I5,I6,I7,I8,I9,I10 DOUBLE PRECISION X(0:NP),Y(0:NP),D0,DNP,D20,D2NP,EPS,BETA,BETAI, * RHO,RHOI,D(0:NP),D2(0:NP),WORK(1:NWORK) C Assign the success value to the error flag. ERRC=0 C Check the size of the work area. IF ( NWORK.LT.COMM+(PART+7)*NP+(N*(N+11))/2+9 ) THEN ERRC=1 RETURN END IF C Compute indices for the splitting of the work array WORK. I1=1 I2=I1+(COMM+PART*NP+NP+1) I3=I2+2*(NP+1) I4=I3+2*(NP+1) I5=I4+(N*(N+1)/2+N) I6=I5+(N+1) I7=I6+(N+1) I8=I7+(N+1) I9=I8+(N+1) I10=I9+NP-1 C DBVSSC is essentially an interfacing routine which relieves the C user of a longer calling sequence. The structure of the method can C be seen in DBVC and in the subroutines called. CALL DBVC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,WORK(I1), * COMM,PART,EPS,KMAX,MAXSTP,WORK(I2),WORK(I3), * WORK(I9),WORK(I10),BETA,BETAI,RHO,RHOI, * D,D2,ERRC,DIAGN) RETURN END C --------------------------------------------------------------------- SUBROUTINE DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C ------------------------------------------------- C Lines 621-754 are comment lines. C The actual code begins at line 760. C ------------------------------------------------- C ABSTRACT: C C DBVSSE is designed to evaluate the interpolating, shape-preserving C spline computed in subroutine DBVSSC. C C C REMARK: C C In these comments variable and array names will be denoted with C capital letters, and with small letters their contents. C C C METHOD: C C Let a spline s:=s(x) of degree n and continuity k (k=1,2) , C interpolating at the knots the point (x(i),y(i)) , i=0,1,...,np , C be previously computed in subroutine DBVSSC. Then, given a set of C tabulation points xtab(i) , i=0,1,...,ntab , DBVSSE computes the C values y0tab(itab):=s(xtab(itab)) and/or C y1tab(itab):=Ds(xtab(itab)) and/or y2tab(itab):=DDs(xtab(itab)) , C using, under user selection, a sequential or binary search scheme. C C The code has the following structure: C C DBVSSE C DBVE C DTRMB C DSQTAB C DLSPIS C DBL C DBL1 C DBL2 C DBNTAB C DBSEAR C DLSPIS C DBL C DBL1 C DBL2 C C C CALLING SEQUENCE: C C CALL DBVSSE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, C * ERRC,D,D2,Y0TAB,Y1TAB,Y2TAB,ERRE,WORK,NWORK) C C C INPUT PARAMETERS: C C X : floating array, of bounds 0:NP, containing the data C abscissas x(i), i=0,1,...,np. Restriction: C x(i).LT.x(i+1), i=0,1,...,np , checked in DBVSSC. C Y : floating array, of bounds 0:NP, containing the data C ordinates y(i), i=0,1,...,np. C NP : integer variable, defining the number of interpolation C points. Restriction: np.GE.2 , checked in DBVSSC. C N : integer variable, containing the degree of s. C Restriction: n.GE.3 , checked in DBVSSC C K : integer variable, containing the class of continuity of s. C Restriction: k.EQ.1 or k.EQ.2 and n.GE.3*k , checked C in DBVSSC. C XTAB : floating array, of bounds 0:NTAB, containing the abscissas C of tabulation points. C Restriction: xtab(i).LE.xtab(i+1), i=0,1,...,ntab-1 . C NTAB : integer variable, defining the number of tabulation points. C Restriction: ntab.GE.0 . C SBOPT : integer variable, containing a control parameter. C If sbopt=1 then the sequential search is used for selecting C the interval of interpolation points in which xtab(i) C falls. If sbopt=2, binary search is used. C Restriction: sbopt.EQ.1 .OR. sbopt.EQ.2 . C Y0OPT : integer variable, containing a control parameter. C If y0opt=1, the spline is evaluated at the points C xtab(i), i=0,1,...,ntab and the results are stored at the C array Y0TAB. C Restriction: y0opt.EQ.0 .OR. y0opt.EQ.1 . C Y1OPT : integer variable, containing a control parameter. C If y1opt=1 the first derivatives of the spline at points C xtab(i) i=0,1,...,ntab , are computed and the results are C stored in the array Y1TAB . C Restriction: y1opt.EQ.0 .OR. y1opt.EQ.1 . C Y2OPT : integer variable, containing a control parameter. C If y2opt=1 the second derivatives of the spline at points C xtab(i), i=0,1,...,ntab are computed and the results are C stored in the array Y2TAB. C Restriction: y2opt.EQ.0 .OR. y2opt.EQ.1 . C ERRC : integer variable, containing the status of the last C execution of subroutine DBVSSC. C D : floating array, of bounds 0:NP, containing the first C derivatives at the knots. C D2 : floating array of bounds 0:NP containing the second C derivatives at the knots. C C C OUTPUT PARAMETERS: C C C Y0TAB : floating array, of bounds 0:NTAB, containing the values of C the spline at the tabulation points xtab(i) , C i=0,1,...,ntab when the option y0opt=1 is activated. C Y1TAB : floating array, of bounds 0:NTAB, containing the values of C the first derivative of the spline at the tabulation points C xtab(i) , i=0,1,...ntab , when the option y1opt=1 is C activated. C Y2TAB : floating array, of bounds 0:NTAB, containing the values of C the second derivative of the spline at the tabulation C points xtab(i) , i=0,1,...,ntab , when the option y2opt=1 C is activated. C ERRE : integer variable, containing an error flag which displays C the status of the code. DBVSSE has only two levels of error C (see DBVSSC for comparison): success and severe error, C which means that some incorrect assignment for input data C have been set. C ERRE=0: success, normal return of the code; C ERRE=1: severe error, the value errc gives a status of C error, which means that the output of DBVSSC is C meaningless. Check the input parameters of DBVSSC. C ERRE=2: severe error, incorrect assignment for some of C the values ntab, sbopt, y0opt, y1opt, y2opt , C nwork; C ERRE=3: severe error, the restriction xtab(i).LT.xtab(i+1) C is not fulfilled for some i when sequential search C is required; C C C OTHER PARAMETERS: C C WORK : floating array, of bounds 1:NKORK, which is used as C a work area to store intermediate results. C The same array can be used to provide workspace for both C the main subroutines DBVSSC and DBVSSE . C NWORK : integer variable containing the size of the work area. C Restriction: nwork .GE. comm+(part+7)*np+(n*(n+11))/2+9 C that is C nwork .GE. 3+(2+7)*np+(n*(n+11))/2+9 C ------------------------------------------------- C End of comments. C ------------------------------------------------- INTEGER COMM,PART PARAMETER (COMM=5, PART=2) INTEGER NP,N,K,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT,ERRC,ERRE,NWORK, * I1,I2,I3,I4,I5,I6,I7,I8,I9,I10 DOUBLE PRECISION X(0:NP),Y(0:NP),XTAB(0:NTAB),Y0TAB(0:NTAB), * Y1TAB(0:NTAB),Y2TAB(0:NTAB),D(0:NP),D2(0:NP), * WORK(1:NWORK) C Assign the success value to the error flag. ERRE=0 C Check the size of the work area. IF ( NWORK.LT.COMM+(PART+7)*NP+(N*(N+11))/2+9 ) THEN ERRE=2 RETURN END IF C Compute indices for the splitting of the work array WORK. I1=1 I2=I1+(COMM+PART*NP+NP+1) I3=I2+2*(NP+1) I4=I3+2*(NP+1) I5=I4+(N*(N+1)/2+N) I6=I5+(N+1) I7=I6+(N+1) I8=I7+(N+1) I9=I8+(N+1) I10=I9+NP-1 C DBVSSE is essentially an interfacing routine which relieves the C user of a longer calling sequence. The structure of the method can C be seen in DBVE and in the subroutines called. CALL DBVE (X,Y,NP,N,K,XTAB,NTAB,SBOPT,Y0OPT,Y1OPT,Y2OPT, * D,D2,ERRC,WORK(I4),WORK(I5),WORK(I6), * WORK(I7),WORK(I8),Y0TAB,Y1TAB,Y2TAB,ERRE) RETURN END SUBROUTINE DALG1 (A1,NP,INFO,COMM,PART,EPS,A2,ERRC,DIAGN) C DALG1 implements the algorithm A1[B(0)] described in subr. DBVSSC. C C The input parameters NP,COMM,PART,EPS and the output parameters C ERRC, DIAGN are described in DBVSSC. The input parameter INFO is C described in DSTINF. C C Items of possible interest are: C C A1: floating array, of bounds 1:2, 0:NP, containing the sequence of C the sets B(i), i=0,1,...,np (see the comments in DBVSSC). C More precisely, B(i) = [a1(1,i),a1(2,i)] . C C A2: floating array, of bounds 1:2, 0:NP, containing the sequence of C the sets A(i), i=0,1,...,np (see the comments in DBVSSC). C More precisely, A(i) = [a2(1,i),a2(2,i)] . INTEGER NP,COMM,PART,ERRC,DIAGN(0:NP-1),ERRC1,I DOUBLE PRECISION A1(1:2,0:NP),INFO(1:COMM+PART*NP+NP+1),EPS, * A2(1:2,0:NP) SAVE FL0 DOUBLE PRECISION FL0 DATA FL0/0.0D00/ ERRC1=0 C Step 1. A2(1,0)=A1(1,0) A2(2,0)=A1(2,0) C Step 2. DO 20 I=1,NP 10 CONTINUE C Step 2.1. CALL DPRJ1 (A2(1,I-1),A2(2,I-1),A1(1,I),A1(2,I),I,INFO, * COMM,PART,NP,A2(1,I),A2(2,I)) C Ignore the constraints in [x(i),x(i+1)] when A(i) is empty. IF (A2(1,I).GT.A2(2,I)+EPS) THEN INFO(COMM+1+(I-1))=FL0 ERRC1=1 DIAGN(I-1)=1 GO TO 10 END IF 20 CONTINUE IF (ERRC1.EQ.1 .AND. ERRC.EQ.9) THEN ERRC=10 ELSE IF (ERRC1.EQ.1) THEN ERRC=6 END IF RETURN END C --------------------------------------------------------------------- SUBROUTINE DALG1D (DSTAR,A1,NP,INFO,COMM,PART,EPS,A2,ERRC1) C DALG1D computes the sequence of sets A(i), i=0,1,...,np, implementing C the algorithm A1[{dstar}], that is with A(0)={dstar} (see the com- C ments in subroutine DBVSSC for details). C C The input parameters NP,COMM,PART,EPS are described in DBVSSC; the C input parameter INFO is described in DSTINF; the input parameters A1 C and A2 are described in subprogram DALG1. C C Item of possible interest is: C C ERRC1 : Integer parameter, containing a control variable which is C then used in subr. DFPSVF C errc1 = 0 - success, normal return of the subprogram; C errc1 = 1 - A(i) is empty for some i. INTEGER NP,COMM,PART,ERRC1,I DOUBLE PRECISION DSTAR,A1(1:2,0:NP),INFO(1:COMM+PART*NP+NP+1), * EPS,A2(1:2,0:NP) ERRC1=0 C Step 1. A2(1,0)=DSTAR A2(2,0)=DSTAR C Step 2. DO 10 I=1,NP CALL DPRJ1 (A2(1,I-1),A2(2,I-1),A1(1,I),A1(2,I),I,INFO, * COMM,PART,NP,A2(1,I),A2(2,I)) IF (A2(1,I).GT.A2(2,I)+EPS) THEN ERRC1=1 RETURN END IF 10 CONTINUE RETURN END C --------------------------------------------------------------------- SUBROUTINE DALG3 (INFO,NP,COMM,PART,OPT,D0,DNP,EPS,KMAX,MAXSTP, * BETA,BETAI,A1,A2,D,ERRC,DIAGN) C DALG3 computes a sequence of slopes ( d(0), d(1), ..., d(np) ) which C can be used to compute a shape-preserving interpolating spline with C or without boundary conditions, as requested by the user. It is an C implementation of the algorithm A3 described in subroutine DBVSSC. C C The input parameters NP,COMM,PART,OPT,EPS,KMAX,MAXSTP,BETA,BETAI,D C and the output parameter ERRC are described in subprogram DBVSSC. C The input parameter INFO is described in subprogram DSTINF. EXTERNAL BETA,BETAI INTEGER NP,COMM,PART,OPT,KMAX,MAXSTP,ERRC,DIAGN(0:NP-1), * I,K,P,Q DOUBLE PRECISION INFO(1:COMM+PART*NP+NP+1),D0,DNP,EPS,BETA, * BETAI,A1(1:2,0:NP),A2(1:2,0:NP),D(0:NP),DSTAR, * P1,P2 LOGICAL DTST INTEGER INK,INSTP PARAMETER (INK=10,INSTP=10) SAVE FL2 DOUBLE PRECISION FL2 DATA FL2/2.0D00/ P=OPT/100 Q=MOD(OPT,100)/10 C If kmax.LE.0 it is reset to ink. IF(KMAX.LE.0) THEN KMAX=INK END IF C If maxstp.LE.0 it is reset to instp. IF(MAXSTP.LE.0) THEN MAXSTP=INSTP END IF C Start step 1: store the sets B(i), i=0,1,...,np , into the array A1. DO 10 I=1,NP CALL DPRJ0 (I,INFO,COMM,PART,NP,A1(1,I-1),A1(2,I-1)) 10 CONTINUE A1(1,NP)=INFO(4) A1(2,NP)=INFO(5) C Reset the first and the last interval if separable boundary condtions C are required IF (Q.EQ.3) THEN A1(1,0)=D0 A1(2,0)=D0 A1(1,NP)=DNP A1(2,NP)=DNP END IF C Start step 2. Call DALG1 to compute the array A2 containing the C sets A(i) , i=0,1,...,np. CALL DALG1 (A1,NP,INFO,COMM,PART,EPS,A2,ERRC,DIAGN) C Start step 3 (steps 3-7 are activated only if boundary conditions are C required). IF (Q.EQ.2) THEN C Compute betai(phi(A(0).INT.beta(A(0))) . CALL DINTRS (MIN(BETA(A2(1,0)),BETA(A2(2,0))), * MAX(BETA(A2(1,0)),BETA(A2(2,0))), * A2(1,NP),A2(2,NP),P1,P2) A1(1,0)=MIN(BETAI(P1),BETAI(P2)) A1(2,0)=MAX(BETAI(P1),BETAI(P2)) C Start step 4. IF (P1.GT.P2+EPS) THEN ERRC=5 RETURN END IF C Start step 5 : initialization K=1 DO 20 I=1,NP A1(1,I)=A2(1,I) A1(2,I)=A2(2,I) 20 CONTINUE C Iteration. The loop is stopped if a convergence test is satisfied C or kmax iterations have already been done. 30 CONTINUE IF (DTST(A1(1,0),A1(2,0),A2(1,0),A2(2,0),EPS).OR. * K.EQ.KMAX) GO TO 50 C Step 5.1 . CALL DALG1 (A1,NP,INFO,COMM,PART,EPS,A2,ERRC,DIAGN) C Step 5.2 . CALL DINTRS (MIN(BETA(A2(1,0)),BETA(A2(2,0))), * MAX(BETA(A2(1,0)),BETA(A2(2,0))), * A2(1,NP),A2(2,NP),P1,P2) A1(1,0)=MIN(BETAI(P1),BETAI(P2)) A1(2,0)=MAX(BETAI(P1),BETAI(P2)) C If gamma(A(0)) is empty for some k the problem does not have any C solution. IF (P1.GT.P2+EPS) THEN ERRC=5 RETURN END IF DO 40 I=1,NP A1(1,I)=A2(1,I) A1(2,I)=A2(2,I) 40 CONTINUE K=K+1 GO TO 30 50 CONTINUE C Start step 7. C Assign to dstar a suitable value IF(P.EQ.1) THEN DSTAR=(A1(1,0)+A1(2,0))/FL2 ELSE DSTAR=INFO(COMM+PART*NP+1) END IF C Check if dstar solves the problem, that is, beta(dstar) belongs to C phi(dstar); if it is not the case, another value for dstar C is looked for. CALL DFPSVF (A1,A2,NP,INFO,COMM,PART,EPS,MAXSTP, * BETA,ERRC,DSTAR) IF(ERRC.NE.0.AND.ERRC.NE.6.AND.ERRC.NE.9 * .AND.ERRC.NE.10) RETURN INFO(COMM+PART*NP+NP+1)=BETA(DSTAR) A2(1,NP)=BETA(DSTAR) A2(2,NP)=BETA(DSTAR) END IF C Start step 8. IF(P.EQ.1) THEN CALL DAL2 (A2,NP,INFO,COMM,PART,D) ELSE CALL DAL2DP (A2,NP,INFO,COMM,PART,D) END IF RETURN END C --------------------------------------------------------------------- SUBROUTINE DAL2(A2,NP,INFO,COMM,PART,D) C DAL2 computes a sequence of slopes (d(0),d(1),...,d(np)) implementing C alg. A2 described in subr. DBVSSC. Each d(i),i=0,1,...,np , is C chosen as the midpoint of the interval of all feasible values . C C The input parameters NP,COMM,PART and the output parameter D are C described in DBVSSC; the input parameter INFO is described in DSTINF. C C Item of possible interest is: C C A2 : floating array, of bounds 1:2, 0:NP; [a2(1,i),a2(2,i)] C is the feasible interval for d(i) . INTEGER NP,COMM,PART,I DOUBLE PRECISION A2(1:2,0:NP),INFO(1:COMM+PART*NP+NP+1), * D(0:NP),P1,P2 SAVE FL1D2 DOUBLE PRECISION FL1D2 DATA FL1D2/0.5D00/ D(NP)=(A2(1,NP)+A2(2,NP))*FL1D2 DO 10 I=NP,1,-1 CALL DPRJ2(A2(1,I-1),A2(2,I-1),D(I),D(I),I,INFO, * COMM,PART,NP,P1,P2) D(I-1)=(P1+P2)*FL1D2 10 CONTINUE RETURN END C --------------------------------------------------------------------- SUBROUTINE DAL2DP (A2,NP,INFO,COMM,PART,D) C DAL2DP links algorithm A2 and a dynamic programming scheme C to select, among the set of all feasible solutions, the sequence C ( d(0),d(1), ..., d(np) ) which is the best 2-norm approximation to C a set of optimal values. More precisely, if (ds(0),ds(1), ...,ds(np)) C is the sequence of optimal values, DAL2DP use the following dynamic C programming relations C C psi(0;d(0)) := (d(0)-ds(0))**2 C psi(j;d(j)) := (d(j)-ds(j))**2 + min(psi(j-1;d(j-1))) C C for describing the objective function C C SUM ((d(j) - ds(j)) ** 2 C j=0,np C C For a complete comprehension of the algorithm see the book \2\ C quoted in the references of subr. DBVSSC C C The input parameters NP,COMM,PART and the output parameter D are C described in subprogram DBVSSC; the input parameter INFO is described C in subprogram DSTINF and the input parameter A2 is described in DAL2. C The constant NSUBD defined below is related to the discretization of C the admissible domain. INTEGER NSUBD PARAMETER (NSUBD=20) INTEGER NP,COMM,PART,IND,I,J,JD0,DMNIND DOUBLE PRECISION A2(1:2,0:NP),INFO(1:COMM+PART*NP+NP+1),D(0:NP), * PSI0(0:NSUBD+1),PSI1(0:NSUBD+1), * PART0(0:NSUBD+1),PART1(0:NSUBD+1),H0,H1,PSI1MN, * P1,P2,D0,D1,DSL SAVE FL0,FLE30 DOUBLE PRECISION FL0,FLE30 DATA FL0,FLE30/0.0D00,1.0D30/ IND=COMM+PART*NP+1 H0=MAX(FL0,(A2(2,0)-A2(1,0))/NSUBD) DO 5 J=0,NSUBD PART0(J)=A2(1,0)+J*H0 5 CONTINUE PART0(NSUBD+1)=DSL(A2(1,0),A2(2,0),INFO(IND)) D(0)=DSL(A2(1,0),A2(2,0),INFO(IND)) DO 10 J=0,NSUBD+1 PSI0(J)=(PART0(J)-D(0))**2 10 CONTINUE DO 40 I=1,NP H1=MAX(FL0,(A2(2,I)-A2(1,I))/NSUBD) DO 15 J=0,NSUBD PART1(J)=A2(1,I)+J*H1 15 CONTINUE PART1(NSUBD+1)=DSL(A2(1,I),A2(2,I),INFO(IND+I)) PSI1MN=FLE30 DO 20 J=0,NSUBD+1 D1=PART1(J) CALL DPRJ2 (A2(1,I-1),A2(2,I-1),D1,D1,I,INFO, * COMM,PART,NP,P1,P2) D0=DSL(P1,P2,D(I-1)) IF (H0.GT.FL0) THEN JD0=DMNIND(D0,PART0) ELSE JD0=0 END IF PSI1(J)=(D1-INFO(IND+I))**2+PSI0(JD0) IF (PSI1(J).LT.PSI1MN) THEN PSI1MN=PSI1(J) D(I)=D1 END IF 20 CONTINUE H0=H1 DO 30 J=0,NSUBD+1 PSI0(J)=PSI1(J) PART0(J)=PART1(J) 30 CONTINUE 40 CONTINUE DO 50 I=NP,1,-1 CALL DPRJ2 (A2(1,I-1),A2(2,I-1),D(I),D(I),I,INFO, * COMM,PART,NP,P1,P2) D(I-1)=DSL(P1,P2,D(I-1)) 50 CONTINUE RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DBL(X,N,L,X0,XN,TB,FLAG,LAUX0) C DBL computes the value assumed by the n-degree Bernstein polynomial C of a function l in the interval (x0,xn) at the point x . C The evaluation is made using a Horner scheme, and the instructions C which do not depend upon x are executed under the control of C the variable FLAG , for avoiding useless computations in C subsequent calls. C The degree n is supposed greater or equal to 3 . C C C INPUT PARAMETERS C C X : floating variable, containing the evaluation point. C N : integer variable, containing the degree of Bernstein C polynomial. C L : floating array, of bounds 0:N , containing the values C of the function l . C X0 : floating variable, containing the left extreme of the C interval. C XN : floating variable, containing the right extreme of the C interval. C TB : floating array, of bounds 0:N , containing the binomial C terms used for computing the Bernstein polynomial. C FLAG : integer variable, containing a control parameter. C In the case flag=0 DBL assumes to perform the first C evaluation of the polynomial, and computes the values C tb(i)*l(i) , i=0,1,...,n . In the case flag=1 DBL C assumes to perform subsequent evaluations, and uses the C values previously computed. C C C OTHER PARAMETERS C C LAUX0 : floating array, of bounds 0:N used as a work area to store C intermediate results. INTEGER N,FLAG,I DOUBLE PRECISION X,L(0:N),X0,XN,TB(0:N),LAUX0(0:N),XNMX,XMX0, * AUX,FL1 SAVE FL1 DATA FL1/1.D00/ IF(FLAG.EQ.0) THEN DO 10 I=0,N LAUX0(I)=TB(I)*L(I) 10 CONTINUE END IF XNMX=XN-X XMX0=X-X0 AUX=FL1 DBL=LAUX0(N) DO 20 I=N-1,0,-1 AUX=XNMX*AUX DBL=LAUX0(I)*AUX+XMX0*DBL 20 CONTINUE DBL=DBL/(XN-X0)**N RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DBL1(X,N,L,X0,XN,TB,FLAG,LAUX1) C DBL1 computes the value assumed by the first derivative of an C n-degree Bernstein polynomial of a function l in the interval C (x0,xn) at the point x . C The evaluation is made using a Horner scheme, and the instructions C which do not depend upon x are executed under the control of C the variable FLAG , for avoiding useless computations in C subsequent calls. C The degree n is supposed greater or equal to 3 . C C INPUT PARAMETERS C C X : floating variable, containing the evaluation point. C N : integer variable, containing the degree of Bernstein C polynomial. C L : floating array, of bounds 0:N , containing the values C of the function l . C X0 : floating variable, containing the left extreme of the C interval. C XN : floating variable, containing the right extreme of the C interval. C TB : floating array, of bounds 0:N-1 , containing the binomial C terms used for computing the Bernstein polynomial. C FLAG : integer variable, containing a control parameter. C In the case flag=0 DBL1 assumes to perform the first C evaluation of the polynomial, and computes the values C tb(i)*(l(i+1)-l(i)) , i=0,1,...,n-1 . In the case flag=1 C DBL1 assumes to perform subsequent evaluations, and uses C the values previously computed. C C C OTHER PARAMETERS C C LAUX1 : floating array, of bounds 0:N-1 used as a work area to store C intermediate results. INTEGER N,FLAG,I DOUBLE PRECISION X,L(0:N),X0,XN,TB(0:N-1),LAUX1(0:N-1),XNMX, * XMX0,AUX,FL1 SAVE FL1 DATA FL1/1.D00/ IF(FLAG.EQ.0) THEN DO 10 I=0,N-1 LAUX1(I)=TB(I)*(L(I+1)-L(I)) 10 CONTINUE END IF XNMX=XN-X XMX0=X-X0 AUX=FL1 DBL1=LAUX1(N-1) DO 20 I=N-2,0,-1 AUX=XNMX*AUX DBL1=LAUX1(I)*AUX+XMX0*DBL1 20 CONTINUE DBL1=N*DBL1/(XN-X0)**N RETURN END C --------------------------------------------------------------------- DOUBLE PRECISION FUNCTION DBL2(X,N,L,X0,XN,TB,FLAG,LAUX2) C DBL2 computes the value assumed by the second derivative of an C n-degree Bernstein polynomial of a function l in the interval C (x0,xn) at the point x . C The evaluation is made using a Horner scheme, and the instructions C which do not depend upon x are executed under the control of C the variable FLAG , for avoiding useless computations in C subsequent calls. C The degree n is supposed greater or equal to 3 . C C INPUT PARAMETERS C C X : floating variable, containing the evaluation point. C N : integer variable, containing the degree of Bernstein C polynomial. C L : floating array, of bounds 0:N , containing the values C of the function l . C X0 : floating variable, containing the left extreme of the C interval. C XN : floating variable, containing the right extreme of the C interval. C TB : floating array, of bounds 0:N-2 , containing the binomial C terms used for computing the Bernstein polynomial. C FLAG : integer variable, containing a control parameter. C In the case flag=0 DBL2 assumes to perform the first C evaluation of the polynomial, and computes the values C tb(i)*(l(i+2)-2*l(i+1)+l(i)) , i=0,1,...,n-2 . C In the case flag=1 DBL2 assumes to perform subsequent C evaluations, and uses the values previously computed. C C C OTHER PARAMETERS C C LAUX2 : floating array, of bounds 0:N-2 used as a work area to store C intermediate results. INTEGER N,FLAG,I DOUBLE PRECISION X,L(0:N),X0,XN,TB(0:N-2),LAUX2(0:N-2),XNMX, * XMX0,AUX,FL1 SAVE FL1 DATA FL1/1.D00/ IF(FLAG.EQ.0) THEN DO 10 I=0,N-2 LAUX2(I)=TB(I)*(L(I+2)-L(I+1)-L(I+1)+L(I)) 10 CONTINUE END IF XNMX=XN-X XMX0=X-X0 AUX=FL1 DBL2=LAUX2(N-2) DO 20 I=N-3,0,-1 AUX=XNMX*AUX DBL2=LAUX2(I)*AUX+XMX0*DBL2 20 CONTINUE DBL2=N*(N-1)*DBL2/(XN-X0)**N RETURN END C --------------------------------------------------------------------- SUBROUTINE DBNTAB (X,Y,NP,XTAB,NTAB,Y0OPT,Y1OPT,Y2OPT,N,K,D,D2, * TB,L,LAUX0,LAUX1,LAUX2,Y0TAB,Y1TAB,Y2TAB) C DBNTAB evaluates the spline and/or its first derivative and/or its C second derivative at the points xtab(j) , j=0,1,...,ntab using C a binary search for finding the interval [x(i),x(i+1)] in which C the tabulation point falls. The input (X,Y,NP,XTAB,NTAB,Y0OPT, C Y1OPT,Y2OPT,N,K,D,D2,TB) and the output (Y0TAB,Y1TAB,Y2TAB) C parameters have been explained in subroutine DBVSSE. For the others C see subroutines DTRMB, DLSPIS. INTEGER NP,NTAB,Y0OPT,Y1OPT,Y2OPT,N,K,IND,J DOUBLE PRECISION X(0:NP),Y(0:NP),XTAB(0:NTAB),D(0:NP),D2(0:NP), * TB(1:N*(N+1)/2+N),L(0:N),LAUX0(0:N),LAUX1(0:N), * LAUX2(0:N),Y0TAB(0:NTAB),Y1TAB(0:NTAB), * Y2TAB(0:NTAB),DBL,DBL1,DBL2 DO 40 J=0,NTAB C Call subprogram DBSEAR to compute the index ind such that C x(ind).LE.xtab(j).LT.x(ind+1) . CALL DBSEAR(X,NP,XTAB(J),IND) C Call subprogram DLSPIS to compute the linear shape-preserving C interpolating spline l at C x(ind)+p*(x(ind+1)-x(ind))/n , p=0,1,...,n . CALL DLSPIS(X,Y,D,D2,NP,N,K,IND,L) IF(Y0OPT.EQ.1) THEN C Evaluate the spline at xtab(j) . Y0TAB(J)=DBL(XTAB(J),N,L,X(IND),X(IND+1), * TB(N*(N+1)/2),0,LAUX0) END IF IF(Y1OPT.EQ.1) THEN C Evaluate the first derivative of the spline at xtab(j) . Y1TAB(J)=DBL1(XTAB(J),N,L,X(IND),X(IND+1), * TB((N-1)*N/2),0,LAUX1) END IF IF(Y2OPT.EQ.1) THEN C Evaluate the second derivative of the spline at xtab(j) . Y2TAB(J)=DBL2(XTAB(J),N,L,X(IND),X(IND+1), * TB((N-2)*(N-1)/2),0,LAUX2) END IF 40 CONTINUE RETURN END C --------------------------------------------------------------------- SUBROUTINE DBSEAR(X,NP,XTAB,IND) C Given an ordered set of points (x(i), i=0,1,...,np) and the C point xtab , DBSEAR finds the index ind such that C C x(ind) .LE. xtab .LT. x(ind+1) C C using a standard binary search. DBSEAR sets ind=0 or ind=np-1 C whenever xtab.LT.x(0) or x(np).LE.xtab . C C C INPUT PARAMETERS C C X : floating array, of bounds 0:NP , containing the set of C ordered points. C XTAB : floating variable, containing the point to be processed. C NP : integer variable, defining the number of points of the C ordered set. C C C OUTPUT PARAMETERS C C IND : integer variable, whose value selects the interval in C which the point xtab falls. INTEGER NP,IND,I1,I2,MED DOUBLE PRECISION X(0:NP),XTAB IF(XTAB.LE.X(0)) THEN IND=0 RETURN END IF IF(XTAB.GE.X(NP)) THEN IND=NP-1 RETURN END IF I1=0 I2=NP 10 CONTINUE IF (.NOT.(I1.NE.I2-1)) GO TO 20 MED=(I1+I2)/2 IF (XTAB.LT.X(MED)) THEN I2=MED ELSE IF (XTAB.GT.X(MED)) THEN I1=MED ELSE IND=MED RETURN END IF GO TO 10 20 CONTINUE IND=I1 RETURN END C --------------------------------------------------------------------- SUBROUTINE DBVC (X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,INFO, * COMM,PART,EPS,KMAX,MAXSTP,A1,A2,DAUX2, * DAUX3,BETA,BETAI,RHO,RHOI,D,D2,ERRC,DIAGN) C DBVC checks input parameters and computes the required spline. C C The input parameters X,Y,NP,N,K,OPT,D0,DNP,D20,D2NP,CONSTR,COMM,PART, C EPS,KMAX,MAXSTP,BETA,BETAI,RHO,RHOI and the output parameters C D,D2,ERRC,DIAGN are described in subroutine DBVSSC. C The other parameters are described in the called subprograms. EXTERNAL BETA,BETAI,RHO,RHOI INTEGER NP,N,K,OPT,CONSTR(0:NP-1),COMM,PART,KMAX,MAXSTP, * ERRC,DIAGN(0:NP-1),P,Q,R,I DOUBLE PRECISION X(0:NP),Y(0:NP),D0,DNP,D20,D2NP, * INFO(1:COMM+PART*NP+NP+1),EPS,A1(1:2,0:NP), * A2(1:2,0:NP),D(0:NP),D2(0:NP),BETA,BETAI, * RHO,RHOI,DAUX2(1:NP-1),DAUX3(0:NP-1) C Check the input parameters NP and OPT. P=OPT/100 Q=MOD(OPT,100)/10 R=MOD(OPT,10) IF ( (NP.LT.2) .OR. (P.LT.1.OR.P.GT.4) .OR. * (Q.LT.1.OR.Q.GT.3) .OR. (R.LT.0.OR.R.GT.4) .OR. * (P.EQ.1.AND.R.EQ.0) ) THEN ERRC=1 RETURN END IF C Check the array CONSTR. IF (MOD(OPT,10).EQ.4) THEN DO 10 I=0,NP-1 IF (CONSTR(I).LT.0 .OR. CONSTR(I).GT.3) THEN ERRC=2 RETURN END IF 10 CONTINUE END IF C Check the input parameters N and K. IF ( (K.LT.1 .OR. K.GT.2) .OR. (N.LT.3*K) ) THEN ERRC=3 RETURN END IF C Check the abscissas of the interpolation points. DO 20 I=0,NP-1 IF(X(I+1).LE.X(I)) THEN ERRC=4 RETURN END IF 20 CONTINUE C Call subprogram DSTINF to set the information array INFO. C Initialize the array DIAGN. DO 30 I=0,NP-1 DIAGN(I)=0 30 CONTINUE CALL DSTINF (OPT,D0,DNP,CONSTR,N,K,X,Y,D,NP,COMM,PART,EPS, * BETA,BETAI,DAUX2,DAUX3,INFO,ERRC,DIAGN) C Call subprogram DALG3 to compute the ar