subroutine banfac ( w, nroww, nrow, nbandl, nbandu, iflag ) !*****************************************************************************80 ! !! banfac() factors a banded matrix without pivoting. ! ! Discussion: ! ! BANFAC returns in W the LU-factorization, without pivoting, of ! the banded matrix A of order NROW with (NBANDL+1+NBANDU) bands ! or diagonals in the work array W. ! ! Gauss elimination without pivoting is used. The routine is ! intended for use with matrices A which do not require row ! interchanges during factorization, especially for the totally ! positive matrices which occur in spline calculations. ! ! The matrix storage mode used is the same one used by LINPACK ! and LAPACK, and results in efficient innermost loops. ! ! Explicitly, A has ! ! NBANDL bands below the diagonal ! 1 main diagonal ! NBANDU bands above the diagonal ! ! and thus, with MIDDLE=NBANDU+1, ! A(I+J,J) is in W(I+MIDDLE,J) for I=-NBANDU,...,NBANDL, J=1,...,NROW. ! ! For example, the interesting entries of a banded matrix ! matrix of order 9, with NBANDL=1, NBANDU=2: ! ! 11 12 13 0 0 0 0 0 0 ! 21 22 23 24 0 0 0 0 0 ! 0 32 33 34 35 0 0 0 0 ! 0 0 43 44 45 46 0 0 0 ! 0 0 0 54 55 56 57 0 0 ! 0 0 0 0 65 66 67 68 0 ! 0 0 0 0 0 76 77 78 79 ! 0 0 0 0 0 0 87 88 89 ! 0 0 0 0 0 0 0 98 99 ! ! would appear in the first 1+1+2=4 rows of W as follows: ! ! 0 0 13 24 35 46 57 68 79 ! 0 12 23 34 45 56 67 78 89 ! 11 22 33 44 55 66 77 88 99 ! 21 32 43 54 65 76 87 98 0 ! ! All other entries of W not identified in this way with an ! entry of A are never referenced. ! ! This routine makes it possible to solve any particular linear system ! A*X=B for X by the call ! ! call banslv ( w, nroww, nrow, nbandl, nbandu, b ) ! ! with the solution X contained in B on return. ! ! If IFLAG=2, then one of NROW-1, NBANDL, NBANDU failed to be nonnegative, ! or else one of the potential pivots was found to be zero ! indicating that A does not have an LU-factorization. This ! implies that A is singular in case it is totally positive. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input/output, real ( kind = rk ) W(NROWW,NROW). ! On input, W contains the "interesting" part of a banded ! matrix A, with the diagonals or bands of A stored in the ! rows of W, while columns of A correspond to columns of W. ! On output, W contains the LU-factorization of A into a unit ! lower triangular matrix L and an upper triangular matrix U ! (both banded) and stored in customary fashion over the ! corresponding entries of A. ! ! Input, integer NROWW, the row dimension of the work array W. ! NROWW must be at least NBANDL+1 + NBANDU. ! ! Input, integer NROW, the number of rows in A. ! ! Input, integer NBANDL, the number of bands of A below ! the main diagonal. ! ! Input, integer NBANDU, the number of bands of A above ! the main diagonal. ! ! Output, integer IFLAG, error flag. ! 1, success. ! 2, failure, the matrix was not factored. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nrow integer nroww real ( kind = rk ) factor integer i integer iflag integer j integer k integer middle integer nbandl integer nbandu real ( kind = rk ) pivot real ( kind = rk ) w(nroww,nrow) iflag = 1 if ( nrow < 1 ) then iflag = 2 return end if ! ! W(MIDDLE,*) contains the main diagonal of A. ! middle = nbandu + 1 if ( nrow == 1 ) then if ( w(middle,nrow) == 0.0D+00 ) then iflag = 2 end if return end if ! ! A is upper triangular. Check that the diagonal is nonzero. ! if ( nbandl <= 0 ) then do i = 1, nrow-1 if ( w(middle,i) == 0.0D+00 ) then iflag = 2 return end if end do if ( w(middle,nrow) == 0.0D+00 ) then iflag = 2 end if return ! ! A is lower triangular. Check that the diagonal is nonzero and ! divide each column by its diagonal. ! else if ( nbandu <= 0 ) then do i = 1, nrow - 1 pivot = w(middle,i) if ( pivot == 0.0D+00 ) then iflag = 2 return end if do j = 1, min ( nbandl, nrow-i ) w(middle+j,i) = w(middle+j,i) / pivot end do end do return end if ! ! A is not just a triangular matrix. ! Construct the LU factorization. ! do i = 1, nrow - 1 ! ! W(MIDDLE,I) is the pivot for the I-th step. ! if ( w(middle,i) == 0.0D+00 ) then iflag = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BANFAC - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot encountered in column ', i stop 1 end if ! ! Divide each entry in column I below the diagonal by PIVOT. ! do j = 1, min ( nbandl, nrow-i ) w(middle+j,i) = w(middle+j,i) / w(middle,i) end do ! ! Subtract A(I,I+K)*(I-th column) from (I+K)-th column (below row I). ! do k = 1, min ( nbandu, nrow-i ) factor = w(middle-k,i+k) do j = 1, min ( nbandl, nrow-i ) w(middle-k+j,i+k) = w(middle-k+j,i+k) - w(middle+j,i) * factor end do end do end do ! ! Check the last diagonal entry. ! if ( w(middle,nrow) == 0.0D+00 ) then iflag = 2 end if return end subroutine banslv ( w, nroww, nrow, nbandl, nbandu, b ) !*****************************************************************************80 ! !! BANSLV solves a banded linear system A * X = B factored by BANFAC. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) W(NROWW,NROW). W contains the banded matrix, ! after it has been factored by BANFAC. ! ! Input, integer NROWW, the row dimension of the work array W. ! NROWW must be at least NBANDL+1 + NBANDU. ! ! Input, integer NROW, the number of rows in A. ! ! Input, integer NBANDL, the number of bands of A below the ! main diagonal. ! ! Input, integer NBANDU, the number of bands of A above the ! main diagonal. ! ! Input/output, real ( kind = rk ) B(NROW). ! On input, B contains the right hand side of the system to be solved. ! On output, B contains the solution, X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nrow integer nroww real ( kind = rk ) b(nrow) integer i integer j integer jmax integer middle integer nbandl integer nbandu real ( kind = rk ) w(nroww,nrow) middle = nbandu + 1 if ( nrow == 1 ) then b(1) = b(1) / w(middle,1) return end if ! ! Forward pass: ! ! For I = 1, 2, ..., NROW-1, subtract RHS(I)*(I-th column of L) ! from the right hand side, below the I-th row. ! if ( 0 < nbandl ) then do i = 1, nrow - 1 jmax = min ( nbandl, nrow-i ) do j = 1, jmax b(i+j) = b(i+j) - b(i) * w(middle+j,i) end do end do end if ! ! Backward pass: ! ! For I=NROW, NROW-1,...,1, divide RHS(I) by ! the I-th diagonal entry of U, then subtract ! RHS(I)*(I-th column of U) from right hand side, above the I-th row. ! do i = nrow, 2, -1 b(i) = b(i) / w(middle,i) do j = 1, min ( nbandu, i - 1 ) b(i-j) = b(i-j) - b(i) * w(middle-j,i) end do end do b(1) = b(1) / w(middle,1) return end subroutine bchfac ( w, nbands, nrow, diag ) !*****************************************************************************80 ! !! BCHFAC constructs a Cholesky factorization of a matrix. ! ! Discussion: ! ! The factorization has the form ! ! C = L * D * L' ! ! with L unit lower triangular and D diagonal, for a given matrix C of ! order NROW, where C is symmetric positive semidefinite and banded, ! having NBANDS diagonals at and below the main diagonal. ! ! Gauss elimination is used, adapted to the symmetry and bandedness of C. ! ! Near-zero pivots are handled in a special way. The diagonal ! element C(N,N) = W(1,N) is saved initially in DIAG(N), all N. ! ! At the N-th elimination step, the current pivot element, W(1,N), ! is compared with its original value, DIAG(N). If, as the result ! of prior elimination steps, this element has been reduced by about ! a word length, that is, if W(1,N) + DIAG(N) <= DIAG(N), then the pivot ! is declared to be zero, and the entire N-th row is declared to ! be linearly dependent on the preceding rows. This has the effect ! of producing X(N) = 0 when solving C * X = B for X, regardless of B. ! ! Justification for this is as follows. In contemplated applications ! of this program, the given equations are the normal equations for ! some least-squares approximation problem, DIAG(N) = C(N,N) gives ! the norm-square of the N-th basis function, and, at this point, ! W(1,N) contains the norm-square of the error in the least-squares ! approximation to the N-th basis function by linear combinations ! of the first N-1. ! ! Having W(1,N)+DIAG(N) <= DIAG(N) signifies that the N-th function ! is linearly dependent to machine accuracy on the first N-1 ! functions, therefore can safely be left out from the basis of ! approximating functions. ! ! The solution of a linear system C * X = B is effected by the ! succession of the following two calls: ! ! call bchfac ( w, nbands, nrow, diag ) ! ! call bchslv ( w, nbands, nrow, b, x ) ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input/output, real ( kind = rk ) W(NBANDS,NROW). ! On input, W contains the NBANDS diagonals in its rows, ! with the main diagonal in row 1. Precisely, W(I,J) ! contains C(I+J-1,J), I=1,...,NBANDS, J=1,...,NROW. ! For example, the interesting entries of a seven diagonal ! symmetric matrix C of order 9 would be stored in W as ! 11 22 33 44 55 66 77 88 99 ! 21 32 43 54 65 76 87 98 * ! 31 42 53 64 75 86 97 * * ! 41 52 63 74 85 96 * * * ! Entries of the array not associated with an ! entry of C are never referenced. ! On output, W contains the Cholesky factorization ! C = L*D*L', with W(1,I) containing 1/D(I,I) and W(I,J) ! containing L(I-1+J,J), I=2,...,NBANDS. ! ! Input, integer NBANDS, indicates the bandwidth of the ! matrix C, that is, C(I,J) = 0 for NBANDS < abs(I-J). ! ! Input, integer NROW, is the order of the matrix C. ! ! Work array, real ( kind = rk ) DIAG(NROW). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbands integer nrow real ( kind = rk ) diag(nrow) integer i integer imax integer j integer jmax integer n real ( kind = rk ) ratio real ( kind = rk ) w(nbands,nrow) if ( nrow <= 1 ) then if ( 0.0D+00 < w(1,1) ) then w(1,1) = 1.0D+00 / w(1,1) end if return end if ! ! Store the diagonal. ! diag(1:nrow) = w(1,1:nrow) ! ! Factorization. ! do n = 1, nrow if ( w(1,n) + diag(n) <= diag(n) ) then w(1:nbands,n) = 0.0D+00 else w(1,n) = 1.0D+00 / w(1,n) imax = min ( nbands - 1, nrow - n ) jmax = imax do i = 1, imax ratio = w(i+1,n) * w(1,n) do j = 1, jmax w(j,n+i) = w(j,n+i) - w(j+i,n) * ratio end do jmax = jmax - 1 w(i+1,n) = ratio end do end if end do return end subroutine bchslv ( w, nbands, nrow, b ) !*****************************************************************************80 ! !! BCHSLV solves a banded symmetric positive definite system. ! ! Discussion: ! ! The system is of the form: ! ! C * X = B ! ! and the Cholesky factorization of C has been constructed ! by BCHFAC. ! ! With the factorization ! ! C = L * D * L' ! ! available, where L is unit lower triangular and D is diagonal, ! the triangular system ! ! L * Y = B ! ! is solved for Y (forward substitution), Y is stored in B, the ! vector D^(-1)*Y is computed and stored in B, then the ! triangular system L'*X = D^(-1)*Y is solved for X ! (back substitution). ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) W(NBANDS,NROW), the Cholesky factorization for C, ! as computed by BCHFAC. ! ! Input, integer NBANDS, the bandwidth of C. ! ! Input, integer NROW, the order of the matrix C. ! ! Input/output, real ( kind = rk ) B(NROW). ! On input, the right hand side. ! On output, the solution. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbands integer nrow real ( kind = rk ) b(nrow) integer j integer n real ( kind = rk ) w(nbands,nrow) if ( nrow <= 1 ) then b(1) = b(1) * w(1,1) return end if ! ! Forward substitution. ! Solve L*Y = B. ! do n = 1, nrow do j = 1, min ( nbands - 1, nrow - n ) b(j+n) = b(j+n) - w(j+1,n) * b(n) end do end do ! ! Back substitution. ! Solve L'*X = D^(-1)*Y. ! do n = nrow, 1, -1 b(n) = b(n) * w(1,n) do j = 1, min ( nbands - 1, nrow - n ) b(n) = b(n) - w(j+1,n) * b(j+n) end do end do return end subroutine bsplpp ( t, bcoef, n, k, scrtch, break, coef, l ) !*****************************************************************************80 ! !! BSPLPP converts from B-spline to piecewise polynomial form. ! ! Discussion: ! ! The B-spline representation of a spline is ! ( T, BCOEF, N, K ), ! while the piecewise polynomial representation is ! ( BREAK, COEF, L, K ). ! ! For each breakpoint interval, the K relevant B-spline coefficients ! of the spline are found and then differenced repeatedly to get the ! B-spline coefficients of all the derivatives of the spline on that ! interval. ! ! The spline and its first K-1 derivatives are then evaluated at the ! left end point of that interval, using BSPLVB repeatedly to obtain ! the values of all B-splines of the appropriate order at that point. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. ! ! Input, real ( kind = rk ) BCOEF(N), the B spline coefficient sequence. ! ! Input, integer N, the number of B spline coefficients. ! ! Input, integer K, the order of the spline. ! ! Work array, real ( kind = rk ) SCRTCH(K,K). ! ! Output, real ( kind = rk ) BREAK(L+1), the piecewise polynomial breakpoint ! sequence. BREAK contains the distinct points in the sequence T(K:N+1) ! ! Output, real ( kind = rk ) COEF(K,N), with COEF(I,J) = (I-1)st derivative ! of the spline at BREAK(J) from the right. ! ! Output, integer L, the number of polynomial pieces which ! make up the spline in the interval ( T(K), T(N+1) ). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer l integer n real ( kind = rk ) bcoef(n) real ( kind = rk ) biatx(k) real ( kind = rk ) break(*) real ( kind = rk ) coef(k,n) real ( kind = rk ) diff integer i integer j integer jp1 integer left integer lsofar real ( kind = rk ) scrtch(k,k) real ( kind = rk ) t(n+k) lsofar = 0 break(1) = t(k) do left = k, n ! ! Find the next nontrivial knot interval. ! if ( t(left+1) == t(left) ) then cycle end if lsofar = lsofar + 1 break(lsofar+1) = t(left+1) if ( k <= 1 ) then coef(1,lsofar) = bcoef(left) cycle end if ! ! Store the K B-spline coefficients relevant to current knot ! interval in SCRTCH(*,1). ! do i = 1, k scrtch(i,1) = bcoef(left-k+i) end do ! ! For J=1,...,K-1, compute the K-J B-spline coefficients relevant to ! the current knot interval for the J-th derivative by differencing ! those for the (J-1)st derivative, and store in SCRTCH(.,J+1). ! do jp1 = 2, k j = jp1 - 1 do i = 1, k - j diff = t(left+i) - t(left+i-(k-j)) if ( 0.0D+00 < diff ) then scrtch(i,jp1) = ( ( scrtch(i+1,j) - scrtch(i,j) ) / diff ) & * real ( k - j, kind = rk ) end if end do end do ! ! For J=0, ..., K-1, find the values at T(left) of the J+1 ! B-splines of order J+1 whose support contains the current ! knot interval from those of order J (in BIATX ), then combine ! with the B-spline coefficients (in SCRTCH(.,K-J) ) found earlier ! to compute the (K-J-1)st derivative at T(LEFT) of the given ! spline. ! call bsplvb ( t, 1, 1, t(left), left, biatx ) coef(k,lsofar) = scrtch(1,k) do jp1 = 2, k call bsplvb ( t, jp1, 2, t(left), left, biatx ) coef(k+1-jp1,lsofar) = dot_product ( biatx(1:jp1), scrtch(1:jp1,k+1-jp1) ) end do end do l = lsofar return end subroutine bsplvb ( t, jhigh, index, x, left, biatx ) !*****************************************************************************80 ! !! BSPLVB evaluates B-splines at a point X with a given knot sequence. ! ! Discusion: ! ! BSPLVB evaluates all possibly nonzero B-splines at X of order ! ! JOUT = MAX ( JHIGH, (J+1)*(INDEX-1) ) ! ! with knot sequence T. ! ! The recurrence relation ! ! X - T(I) T(I+J+1) - X ! B(I,J+1)(X) = ----------- * B(I,J)(X) + --------------- * B(I+1,J)(X) ! T(I+J)-T(I) T(I+J+1)-T(I+1) ! ! is used to generate B(LEFT-J:LEFT,J+1)(X) from B(LEFT-J+1:LEFT,J)(X) ! storing the new values in BIATX over the old. ! ! The facts that ! ! B(I,1)(X) = 1 if T(I) <= X < T(I+1) ! ! and that ! ! B(I,J)(X) = 0 unless T(I) <= X < T(I+J) ! ! are used. ! ! The particular organization of the calculations follows ! algorithm 8 in chapter X of the text. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(LEFT+JOUT), the knot sequence. T is assumed to ! be nondecreasing, and also, T(LEFT) must be strictly less than ! T(LEFT+1). ! ! Input, integer JHIGH, INDEX, determine the order ! JOUT = max ( JHIGH, (J+1)*(INDEX-1) ) ! of the B-splines whose values at X are to be returned. ! INDEX is used to avoid recalculations when several ! columns of the triangular array of B-spline values are ! needed, for example, in BVALUE or in BSPLVD. ! If INDEX = 1, the calculation starts from scratch and the entire ! triangular array of B-spline values of orders ! 1, 2, ...,JHIGH is generated order by order, that is, ! column by column. ! If INDEX = 2, only the B-spline values of order J+1, J+2, ..., JOUT ! are generated, the assumption being that BIATX, J, ! DELTAL, DELTAR are, on entry, as they were on exit ! at the previous call. In particular, if JHIGH = 0, ! then JOUT = J+1, that is, just the next column of B-spline ! values is generated. ! Warning: the restriction JOUT <= JMAX (= 20) is ! imposed arbitrarily by the dimension statement for DELTAL ! and DELTAR, but is nowhere checked for. ! ! Input, real ( kind = rk ) X, the point at which the B-splines ! are to be evaluated. ! ! Input, integer LEFT, an integer chosen so that ! T(LEFT) <= X <= T(LEFT+1). ! ! Output, real ( kind = rk ) BIATX(JOUT), with BIATX(I) containing the ! value at X of the polynomial of order JOUT which agrees ! with the B-spline B(LEFT-JOUT+I,JOUT,T) on the interval ! (T(LEFT),T(LEFT+1)). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: jmax = 20 integer jhigh real ( kind = rk ) biatx(jhigh) real ( kind = rk ), save, dimension ( jmax ) :: deltal real ( kind = rk ), save, dimension ( jmax ) :: deltar integer i integer index integer, save :: j = 1 integer left real ( kind = rk ) saved real ( kind = rk ) t(left+jhigh) real ( kind = rk ) term real ( kind = rk ) x if ( index == 1 ) then j = 1 biatx(1) = 1.0D+00 if ( jhigh <= j ) then return end if end if if ( t(left+1) <= t(left) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BSPLVB - Fatal error!' write ( *, '(a)' ) ' It is required that T(LEFT) < T(LEFT+1).' write ( *, '(a,i8)' ) ' But LEFT = ', left write ( *, '(a,g14.6)' ) ' T(LEFT) = ', t(left) write ( *, '(a,g14.6)' ) ' T(LEFT+1) = ', t(left+1) stop 1 end if do deltar(j) = t(left+j) - x deltal(j) = x - t(left+1-j) saved = 0.0D+00 do i = 1, j term = biatx(i) / ( deltar(i) + deltal(j+1-i) ) biatx(i) = saved + deltar(i) * term saved = deltal(j+1-i) * term end do biatx(j+1) = saved j = j + 1 if ( jhigh <= j ) then exit end if end do return end subroutine bsplvd ( t, k, x, left, a, dbiatx, nderiv ) !*****************************************************************************80 ! !! BSPLVD calculates the nonvanishing B-splines and derivatives at X. ! ! Discussion: ! ! Values at X of all the relevant B-splines of order K:K+1-NDERIV ! are generated via BSPLVB and stored temporarily in DBIATX. ! ! Then the B-spline coefficients of the required derivatives ! of the B-splines of interest are generated by differencing, ! each from the preceding one of lower order, and combined with ! the values of B-splines of corresponding order in DBIATX ! to produce the desired values. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(LEFT+K), the knot sequence. It is assumed that ! T(LEFT) < T(LEFT+1). Also, the output is correct only if ! T(LEFT) <= X <= T(LEFT+1). ! ! Input, integer K, the order of the B-splines to be evaluated. ! ! Input, real ( kind = rk ) X, the point at which these values are sought. ! ! Input, integer LEFT, indicates the left endpoint of the ! interval of interest. The K B-splines whose support contains the interval ! ( T(LEFT), T(LEFT+1) ) are to be considered. ! ! Workspace, real ( kind = rk ) A(K,K). ! ! Output, real ( kind = rk ) DBIATX(K,NDERIV). DBIATX(I,M) contains ! the value of the (M-1)st derivative of the (LEFT-K+I)-th B-spline ! of order K for knot sequence T, I=M,...,K, M=1,...,NDERIV. ! ! Input, integer NDERIV, indicates that values of B-splines and ! their derivatives up to but not including the NDERIV-th are asked for. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer left integer nderiv real ( kind = rk ) a(k,k) real ( kind = rk ) dbiatx(k,nderiv) real ( kind = rk ) factor real ( kind = rk ) fkp1mm integer i integer ideriv integer il integer j integer jlow integer jp1mid integer ldummy integer m integer mhigh real ( kind = rk ) t(left+k) real ( kind = rk ) x mhigh = max ( min ( nderiv, k ), 1 ) ! ! MHIGH is usually equal to NDERIV. ! call bsplvb ( t, k+1-mhigh, 1, x, left, dbiatx ) if ( mhigh == 1 ) then return end if ! ! The first column of DBIATX always contains the B-spline values ! for the current order. These are stored in column K+1-current ! order before BSPLVB is called to put values for the next ! higher order on top of it. ! ideriv = mhigh do m = 2, mhigh jp1mid = 1 do j = ideriv, k dbiatx(j,ideriv) = dbiatx(jp1mid,1) jp1mid = jp1mid + 1 end do ideriv = ideriv - 1 call bsplvb ( t, k+1-ideriv, 2, x, left, dbiatx ) end do ! ! At this point, B(LEFT-K+I, K+1-J)(X) is in DBIATX(I,J) for ! I=J,...,K and J=1,...,MHIGH ('=' NDERIV). ! ! In particular, the first column of DBIATX is already in final form. ! ! To obtain corresponding derivatives of B-splines in subsequent columns, ! generate their B-representation by differencing, then evaluate at X. ! jlow = 1 do i = 1, k a(jlow:k,i) = 0.0D+00 jlow = i a(i,i) = 1.0D+00 end do ! ! At this point, A(.,J) contains the B-coefficients for the J-th of the ! K B-splines of interest here. ! do m = 2, mhigh fkp1mm = real ( k + 1 - m, kind = rk ) il = left i = k ! ! For J = 1,...,K, construct B-coefficients of (M-1)st derivative of ! B-splines from those for preceding derivative by differencing ! and store again in A(.,J). The fact that A(I,J) = 0 for ! I < J is used. ! do ldummy = 1, k + 1 - m factor = fkp1mm / ( t(il+k+1-m) - t(il) ) ! ! The assumption that T(LEFT) < T(LEFT+1) makes denominator ! in FACTOR nonzero. ! a(i,1:i) = ( a(i,1:i) - a(i-1,1:i) ) * factor il = il - 1 i = i - 1 end do ! ! For I = 1,...,K, combine B-coefficients A(.,I) with B-spline values ! stored in DBIATX(.,M) to get value of (M-1)st derivative of ! I-th B-spline (of interest here) at X, and store in DBIATX(I,M). ! ! Storage of this value over the value of a B-spline ! of order M there is safe since the remaining B-spline derivatives ! of the same order do not use this value due to the fact ! that A(J,I) = 0 for J < I. ! do i = 1, k jlow = max ( i, m ) dbiatx(i,m) = dot_product ( a(jlow:k,i), dbiatx(jlow:k,m) ) end do end do return end subroutine bspp2d ( t, bcoef, n, k, m, scrtch, break, coef, l ) !*****************************************************************************80 ! !! BSPP2D converts from B-spline to piecewise polynomial representation. ! ! Discussion: ! ! The B-spline representation ! ! T, BCOEF(.,J), N, K ! ! is converted to its piecewise polynomial representation ! ! BREAK, COEF(J,.,.), L, K, J=1, ..., M. ! ! This is an extended version of BSPLPP for use with tensor products. ! ! For each breakpoint interval, the K relevant B-spline ! coefficients of the spline are found and then differenced ! repeatedly to get the B-spline coefficients of all the ! derivatives of the spline on that interval. ! ! The spline and its first K-1 derivatives are then evaluated ! at the left endpoint of that interval, using BSPLVB ! repeatedly to obtain the values of all B-splines of the ! appropriate order at that point. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. ! ! Input, real ( kind = rk ) BCOEF(N,M). For each J, B(*,J) is the ! B-spline coefficient sequence, of length N. ! ! Input, integer N, the length of BCOEF. ! ! Input, integer K, the order of the spline. ! ! Input, integer M, the number of data sets. ! ! Work array, real ( kind = rk ) SCRTCH(K,K,M). ! ! Output, real ( kind = rk ) BREAK(L+1), the breakpoint sequence ! containing the distinct points in the sequence T(K),...,T(N+1) ! ! Output, real ( kind = rk ) COEF(M,K,N), with COEF(MM,I,J) = the (I-1)st ! derivative of the MM-th spline at BREAK(J) from the right, MM=1, ..., M. ! ! Output, integer L, the number of polynomial pieces which ! make up the spline in the interval (T(K), T(N+1)). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer m integer n real ( kind = rk ) bcoef(n,m) real ( kind = rk ) biatx(k) real ( kind = rk ) break(*) real ( kind = rk ) coef(m,k,*) real ( kind = rk ) diff real ( kind = rk ) fkmj integer i integer j integer jp1 integer kmj integer l integer left integer lsofar integer mm real ( kind = rk ) scrtch(k,k,m) real ( kind = rk ) sum1 real ( kind = rk ) t(n+k) lsofar = 0 break(1) = t(k) do left = k, n ! ! Find the next nontrivial knot interval. ! if ( t(left+1) == t(left) ) then cycle end if lsofar = lsofar + 1 break(lsofar+1) = t(left+1) if ( k <= 1 ) then coef(1:m,1,lsofar) = bcoef(left,1:m) cycle end if ! ! Store the K B-spline coefficients relevant to current knot interval ! in SCRTCH(.,1). ! do i = 1, k scrtch(i,1,1:m) = bcoef(left-k+i,1:m) end do ! ! For J = 1,...,K-1, compute the ( K - J ) B-spline coefficients relevant to ! current knot interval for the J-th derivative by differencing ! those for the (J-1)st derivative, and store in SCRTCH(.,J+1). ! do jp1 = 2, k j = jp1 - 1 kmj = k - j fkmj = real ( k - j, kind = rk ) do i = 1, k - j diff = ( t(left+i) - t(left+i-kmj) ) / fkmj if ( 0.0D+00 < diff ) then scrtch(i,jp1,1:m) = ( scrtch(i+1,j,1:m) - scrtch(i,j,1:m) ) / diff end if end do end do ! ! For J = 0, ..., K-1, find the values at T(LEFT) of the J+1 ! B-splines of order J+1 whose support contains the current ! knot interval from those of order J (in BIATX ), then combine ! with the B-spline coefficients (in SCRTCH(.,K-J) ) found earlier ! to compute the (K-J-1)st derivative at T(LEFT) of the given spline. ! call bsplvb ( t, 1, 1, t(left), left, biatx ) coef(1:m,k,lsofar) = scrtch(1,k,1:m) do jp1 = 2, k call bsplvb ( t, jp1, 2, t(left), left, biatx ) kmj = k + 1 - jp1 do mm = 1, m sum1 = 0.0D+00 do i = 1, jp1 sum1 = sum1 + biatx(i) * scrtch(i,kmj,mm) end do coef(mm,kmj,lsofar) = sum1 end do end do end do l = lsofar return end function bvalue ( t, bcoef, n, k, x, jderiv ) !*****************************************************************************80 ! !! BVALUE evaluates a derivative of a spline from its B-spline representation. ! ! Discussion: ! ! The spline is taken to be continuous from the right. ! ! The nontrivial knot interval (T(I),T(I+1)) containing X is ! located with the aid of INTERV. The K B-spline coefficients ! of F relevant for this interval are then obtained from BCOEF, ! or are taken to be zero if not explicitly available, and are ! then differenced JDERIV times to obtain the B-spline ! coefficients of (D^JDERIV)F relevant for that interval. ! ! Precisely, with J = JDERIV, we have from X.(12) of the text that: ! ! (D^J)F = sum ( BCOEF(.,J)*B(.,K-J,T) ) ! ! where ! / BCOEF(.), if J == 0 ! / ! BCOEF(.,J) = / BCOEF(.,J-1) - BCOEF(.-1,J-1) ! / -----------------------------, if 0 < J ! / (T(.+K-J) - T(.))/(K-J) ! ! Then, we use repeatedly the fact that ! ! sum ( A(.) * B(.,M,T)(X) ) = sum ( A(.,X) * B(.,M-1,T)(X) ) ! ! with ! (X - T(.))*A(.) + (T(.+M-1) - X)*A(.-1) ! A(.,X) = --------------------------------------- ! (X - T(.)) + (T(.+M-1) - X) ! ! to write (D^J)F(X) eventually as a linear combination of ! B-splines of order 1, and the coefficient for B(I,1,T)(X) ! must then be the desired number (D^J)F(X). ! See Chapter X, (17)-(19) of text. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. T is assumed ! to be nondecreasing. ! ! Input, real ( kind = rk ) BCOEF(N), B-spline coefficient sequence. ! ! Input, integer N, the length of BCOEF. ! ! Input, integer K, the order of the spline. ! ! Input, real ( kind = rk ) X, the point at which to evaluate. ! ! Input, integer JDERIV, the order of the derivative to ! be evaluated. JDERIV is assumed to be zero or positive. ! ! Output, real ( kind = rk ) BVALUE, the value of the (JDERIV)-th ! derivative of the spline at X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer n real ( kind = rk ) aj(k) real ( kind = rk ) bcoef(n) real ( kind = rk ) bvalue real ( kind = rk ) dl(k) real ( kind = rk ) dr(k) integer i integer ilo integer j integer jc integer jcmax integer jcmin integer jderiv integer jj integer mflag real ( kind = rk ) t(n+k) real ( kind = rk ) x bvalue = 0.0D+00 if ( k <= jderiv ) then return end if ! ! Find I so that 1 <= I < N+K and T(I) < T(I+1) and T(I) <= X < T(I+1). ! ! If no such I can be found, X lies outside the support of the ! spline F and BVALUE = 0. The asymmetry in this choice of I makes F ! right continuous, except at T(N+K) where it is leftcontinuous. ! call interv ( t, n+k, x, i, mflag ) if ( mflag /= 0 ) then return end if ! ! If K = 1 (and JDERIV = 0), BVALUE = BCOEF(I). ! if ( k <= 1 ) then bvalue = bcoef(i) return end if ! ! Store the K B-spline coefficients relevant for the knot interval ! ( T(I),T(I+1) ) in AJ(1),...,AJ(K) and compute DL(J) = X - T(I+1-J), ! DR(J) = T(I+J)-X, J=1,...,K-1. Set any of the AJ not obtainable ! from input to zero. ! ! Set any T's not obtainable equal to T(1) or to T(N+K) appropriately. ! jcmin = 1 if ( k <= i ) then do j = 1, k-1 dl(j) = x - t(i+1-j) end do else jcmin = 1 - ( i - k ) do j = 1, i dl(j) = x - t(i+1-j) end do do j = i, k-1 aj(k-j) = 0.0D+00 dl(j) = dl(i) end do end if jcmax = k if ( n < i ) then jcmax = k + n - i do j = 1, k + n - i dr(j) = t(i+j) - x end do do j = k+n-i, k-1 aj(j+1) = 0.0D+00 dr(j) = dr(k+n-i) end do else do j = 1, k-1 dr(j) = t(i+j) - x end do end if do jc = jcmin, jcmax aj(jc) = bcoef(i-k+jc) end do ! ! Difference the coefficients JDERIV times. ! do j = 1, jderiv ilo = k - j do jj = 1, k - j aj(jj) = ( ( aj(jj+1) - aj(jj) ) / ( dl(ilo) + dr(jj) ) ) & * real ( k - j, kind = rk ) ilo = ilo - 1 end do end do ! ! Compute value at X in (T(I),T(I+1)) of JDERIV-th derivative, ! given its relevant B-spline coefficients in AJ(1),...,AJ(K-JDERIV). ! do j = jderiv+1, k-1 ilo = k-j do jj = 1, k-j aj(jj) = ( aj(jj+1) * dl(ilo) + aj(jj) * dr(jj) ) & / ( dl(ilo) + dr(jj) ) ilo = ilo - 1 end do end do bvalue = aj(1) return end subroutine chol1d ( p, v, qty, npoint, ncol, u, qu ) !*****************************************************************************80 ! !! CHOL1D sets up and solves linear systems needed by SMOOTH. ! ! Discussion: ! ! This routine constructs the upper three diagonals of ! ! V(I,J), I = 2 to NPOINT-1, J=1,3, ! ! of the matrix ! ! 6 * (1-P) * Q' * (D^2) * Q + P * R. ! ! It then computes its L*L' decomposition and stores it also ! in V, then applies forward and back substitution to the right hand side ! ! Q'*Y ! ! in QTY to obtain the solution in U. ! ! Modified: ! ! 16 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) P, the smoothing parameter that defines ! the linear system. ! ! Input/output, real ( kind = rk ) V(NPOINT,7), contains data used ! to define the linear system, some of which is determined by ! routine SETUPQ. ! ! Input, real ( kind = rk ) QTY(NPOINT), the value of Q' * Y. ! ! Input, integer NPOINT, the number of equations. ! ! Input, integer NCOL, an unused parameter, which may be ! set to 1. ! ! Output, real ( kind = rk ) U(NPOINT), the solution. ! ! Output, real ( kind = rk ) QU(NPOINT), the value of Q * U. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer npoint integer i integer ncol real ( kind = rk ) p real ( kind = rk ) qty(npoint) real ( kind = rk ) qu(npoint) real ( kind = rk ) u(npoint) real ( kind = rk ) v(npoint,7) real ( kind = rk ) prev real ( kind = rk ) ratio real ( kind = rk ) six1mp real ( kind = rk ) twop ! ! Construct 6*(1-P)*Q'*(D^2)*Q + P*R. ! six1mp = 6.0D+00 * ( 1.0D+00 - p ) twop = 2.0D+00 * p v(2:npoint-1,1) = six1mp * v(2:npoint-1,5) & + twop * ( v(1:npoint-2,4) + v(2:npoint-1,4) ) v(2:npoint-1,2) = six1mp * v(2:npoint-1,6) + p * v(2:npoint-1,4) v(2:npoint-1,3) = six1mp * v(2:npoint-1,7) if ( npoint < 4 ) then u(1) = 0.0D+00 u(2) = qty(2) / v(2,1) u(3) = 0.0D+00 ! ! Factorization. ! else do i = 2, npoint-2 ratio = v(i,2) / v(i,1) v(i+1,1) = v(i+1,1) - ratio * v(i,2) v(i+1,2) = v(i+1,2) - ratio * v(i,3) v(i,2) = ratio ratio = v(i,3) / v(i,1) v(i+2,1) = v(i+2,1) - ratio * v(i,3) v(i,3) = ratio end do ! ! Forward substitution ! u(1) = 0.0D+00 v(1,3) = 0.0D+00 u(2) = qty(2) do i = 2, npoint-2 u(i+1) = qty(i+1) - v(i,2) * u(i) - v(i-1,3) * u(i-1) end do ! ! Back substitution. ! u(npoint) = 0.0D+00 u(npoint-1) = u(npoint-1) / v(npoint-1,1) do i = npoint-2, 2, -1 u(i) = u(i) / v(i,1) - u(i+1) * v(i,2) - u(i+2) * v(i,3) end do end if ! ! Construct Q * U. ! prev = 0.0D+00 do i = 2, npoint qu(i) = ( u(i) - u(i-1) ) / v(i-1,4) qu(i-1) = qu(i) - prev prev = qu(i) end do qu(npoint) = -qu(npoint) return end subroutine colloc ( aleft, aright, lbegin, iorder, ntimes, addbrk, relerr ) !*****************************************************************************80 ! !! COLLOC solves an ordinary differential equation by collocation. ! ! Method: ! ! The M-th order ordinary differential equation with M side ! conditions, to be specified in subroutine DIFEQU, is solved ! approximately by collocation. ! ! The approximation F to the solution G is piecewise polynomial of order ! K+M with L pieces and M-1 continuous derivatives. F is determined by ! the requirement that it satisfy the differential equation at K points ! per interval (to be specified in COLPNT ) and the M side conditions. ! ! This usually nonlinear system of equations for F is solved by ! Newton's method. the resulting linear system for the B-coefficients of an ! iterate is constructed appropriately in EQBLOK and then solved ! in SLVBLK, a program designed to solve almost block ! diagonal linear systems efficiently. ! ! There is an opportunity to attempt improvement of the breakpoint ! sequence, both in number and location, through the use of NEWNOT. ! ! Printed output consists of the piecewise polynomial representation ! of the approximate solution, and of the error at selected points. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) ALEFT, ARIGHT, the endpoints of the interval. ! ! Input, integer LBEGIN, the initial number of polynomial ! pieces in the approximation. A uniform breakpoint sequence will be chosen. ! ! Input, integer IORDER, the order of the polynomial pieces ! to be used in the approximation ! ! Input, integer NTIMES, the number of passes to be made ! through NEWNOT. ! ! Input, real ( kind = rk ) ADDBRK, the number, possibly fractional, of ! breaks to be added per pass through NEWNOT. For instance, if ! ADDBRK = 0.33334, then a breakpoint will be added at every third pass ! through NEWNOT. ! ! Input, real ( kind = rk ) RELERR, a tolerance. Newton iteration is ! stopped if the difference between the B-coefficients of two successive ! iterates is no more than RELERR*(absolute largest B-coefficient). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: npiece = 100 integer, parameter :: ndim = 200 integer, parameter :: ncoef = 2000 integer, parameter :: lenblk = 2000 real ( kind = rk ) a(ndim) real ( kind = rk ) addbrk real ( kind = rk ) aleft real ( kind = rk ) amax real ( kind = rk ) aright real ( kind = rk ) asave(ndim) real ( kind = rk ) b(ndim) real ( kind = rk ) bloks(lenblk) real ( kind = rk ) break real ( kind = rk ) coef real ( kind = rk ) dx real ( kind = rk ) err integer i integer iflag integer ii integer integs(3,npiece) integer iorder integer iside integer itemps(ndim) integer iter integer itermx integer k integer kpm integer l integer lbegin integer lnew integer m integer n integer nbloks integer nt integer ntimes real ( kind = rk ) relerr real ( kind = rk ) rho real ( kind = rk ) t(ndim) real ( kind = rk ) templ(lenblk) real ( kind = rk ) temps(ndim) real ( kind = rk ) xside equivalence ( bloks, templ ) save / approx / save / other / save / side / common / approx / break(npiece), coef(ncoef), l, kpm common / other / itermx, k, rho(19) common / side / m, iside, xside(10) kpm = iorder if ( ncoef < lbegin * kpm ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLLOC - Fatal error!' write ( *, '(a)' ) ' The assigned dimension for COEF is too small.' stop 1 end if ! ! Set the various parameters concerning the particular differential ! equation, including a first approximation in case the differential ! equation is to be solved by iteration ( 0 < ITERMX ). ! call difequ ( 1, temps(1), temps ) ! ! Obtain the K collocation points for the standard interval. ! k = kpm - m call colpnt ( k, rho ) ! ! The following five statements could be replaced by a read in ! order to obtain a nonuniform spacing of the breakpoints. ! dx = ( aright - aleft ) / real ( lbegin, kind = rk ) temps(1) = aleft do i = 2, lbegin temps(i) = temps(i-1) + dx end do temps(lbegin+1) = aright ! ! Generate the required knots T(1:N+KPM). ! call knots ( temps, lbegin, kpm, m, t, n ) nt = 1 ! ! Generate the almost block diagonal coefficient matrix BLOKS and ! right hand side B from collocation equations and side conditions. ! ! Then solve via SLVBLK, obtaining the B-representation of the ! approximation in T, A, N, KPM. ! do call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b ) call slvblk ( bloks, integs, nbloks, b, itemps, a, iflag ) ! ! Save B-spline coefficients of current approximation in ASAVE, then ! get new approximation and compare with old. ! ! If coefficients are more than RELERR apart (relatively) or if number ! of iterations is less than ITERMX, continue iterating. ! do iter = 1, itermx call bsplpp ( t, a, n, kpm, templ, break, coef, l ) asave(1:n) = a(1:n) call eqblok ( t, n, kpm, temps, a, bloks, lenblk, integs, nbloks, b ) call slvblk ( bloks, integs, nbloks, b, itemps, a, iflag ) amax = maxval ( abs ( a(1:n) ) ) err = maxval ( abs ( a(1:n) - asave(1:n) ) ) if ( err <= relerr * amax ) then exit end if end do ! ! Iteration (if any) completed. Print out approximation based on current ! breakpoint sequence, then try to improve the sequence. ! write ( *, '(a)' ) ' ' write ( *,'(a,i3,a,i3,a)' ) & ' Approximation from a space of splines of order ', kpm, & ' on ', l, ' intervals' write ( *, '(a,i4)' ) ' of dimension ', n write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Breakpoints:' write ( *, '(a)' ) ' ' write ( *, '(5g14.6)' ) break(2:l) if ( 0 < itermx ) then write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Results on interation ', iter end if call bsplpp ( t, a, n, kpm, templ, break, coef, l ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) & ' The piecewise polynomial representation of the approximation:' write ( *, '(a)' ) ' ' do i = 1, l ii = ( i - 1 ) * kpm write ( *, '(f9.3,2x,e12.4,10e11.3)' ) break(i), coef(ii+1:ii+kpm) end do ! ! The following call is provided here for possible further analysis ! of the approximation specific to the problem being solved. ! It is, of course, easily omitted. ! call difequ ( 4, temps(1), temps ) if ( ntimes < nt ) then exit end if ! ! From the piecewise polynomial representation of the current approximation, ! obtain in NEWNOT a new, and possibly better, sequence of breakpoints, ! adding, on average, ADDBRK breakpoints per pass through NEWNOT. ! lnew = lbegin + int ( real ( nt, kind = rk ) * addbrk ) if ( ncoef < lnew * kpm ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLLOC - Fatal error!' write ( *, '(a)' ) ' The assigned dimension for COEF is too small.' stop 1 end if call newnot ( break, coef, l, kpm, temps, lnew, templ ) call knots ( temps, lnew, kpm, m, t, n ) nt = nt + 1 end do return end subroutine colpnt ( k, rho ) !*****************************************************************************80 ! !! COLPNT supplies collocation points. ! ! Discussion: ! ! The collocation points are for the standard interval (-1,1) as the ! zeros of the Legendre polynomial of degree K, provided K <= 8. ! ! Otherwise, uniformly spaced points are given. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, integer K, the number of collocation points desired. ! ! Output, real ( kind = rk ) RHO(K), the collocation points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer j real ( kind = rk ) rho(k) if ( k == 1 ) then rho(1) = 0.0D+00 else if ( k == 2 ) then rho(1) = -0.577350269189626D+00 rho(2) = 0.577350269189626D+00 else if ( k == 3 ) then rho(1) = -0.774596669241483D+00 rho(2) = 0.0D+00 rho(3) = 0.774596669241483D+00 else if ( k == 4 ) then rho(1) = -0.861136311594053D+00 rho(2) = -0.339981043584856D+00 rho(3) = 0.339981043584856D+00 rho(4) = 0.861136311594053D+00 else if ( k == 5 ) then rho(1) = -0.906179845938664D+00 rho(2) = -0.538469310105683D+00 rho(3) = 0.0D+00 rho(4) = 0.538469310105683D+00 rho(5) = 0.906179845938664D+00 else if ( k == 6 ) then rho(1) = -0.932469514203152D+00 rho(2) = -0.661209386466265D+00 rho(3) = -0.238619186083197D+00 rho(4) = 0.238619186083197D+00 rho(5) = 0.661209386466265D+00 rho(6) = 0.932469514203152D+00 else if ( k == 7 ) then rho(1) = -0.949107912342759D+00 rho(2) = -0.741531185599394D+00 rho(3) = -0.405845151377397D+00 rho(4) = 0.0D+00 rho(5) = 0.405845151377397D+00 rho(6) = 0.741531185599394D+00 rho(7) = 0.949107912342759D+00 else if ( k == 8 ) then rho(1) = -0.960289856497536D+00 rho(2) = -0.796666477413627D+00 rho(3) = -0.525532409916329D+00 rho(4) = -0.183434642495650D+00 rho(5) = 0.183434642495650D+00 rho(6) = 0.525532409916329D+00 rho(7) = 0.796666477413627D+00 rho(8) = 0.960289856497536D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' )'COLPNT - Warning!' write ( *, '(a)' )' Equispaced collocation points will be used,' write ( *, '(a,i8)' ) ' because K = ', k do j = 1, k rho(j) = ( real ( k - j, kind = rk ) * ( -1.0D+00 ) & + real ( j - 1, kind = rk ) * ( +1.0D+00 ) ) & / real ( k - 1, kind = rk ) end do end if return end subroutine cubspl ( tau, c, n, ibcbeg, ibcend ) !*****************************************************************************80 ! !! CUBSPL defines an interpolatory cubic spline. ! ! Discussion: ! ! A tridiagonal linear system for the unknown slopes S(I) of ! F at TAU(I), I=1,..., N, is generated and then solved by Gauss ! elimination, with S(I) ending up in C(2,I), for all I. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) TAU(N), the abscissas or X values of ! the data points. The entries of TAU are assumed to be ! strictly increasing. ! ! Input, integer N, the number of data points. N is ! assumed to be at least 2. ! ! Input/output, real ( kind = rk ) C(4,N). ! On input, if IBCBEG or IBCBEG is 1 or 2, then C(2,1) ! or C(2,N) should have been set to the desired derivative ! values, as described further under IBCBEG and IBCEND. ! On output, C contains the polynomial coefficients of ! the cubic interpolating spline with interior knots ! TAU(2) through TAU(N-1). ! In the interval interval (TAU(I), TAU(I+1)), the spline ! F is given by ! F(X) = ! C(1,I) + ! C(2,I) * H + ! C(3,I) * H^2 / 2 + ! C(4,I) * H^3 / 6. ! where H=X-TAU(I). The routine PPVALU may be used to ! evaluate F or its derivatives from TAU, C, L=N-1, ! and K=4. ! ! Input, integer IBCBEG, IBCEND, boundary condition indicators. ! IBCBEG = 0 means no boundary condition at TAU(1) is given. ! In this case, the "not-a-knot condition" is used. That ! is, the jump in the third derivative across TAU(2) is ! forced to zero. Thus the first and the second cubic ! polynomial pieces are made to coincide. ! IBCBEG = 1 means the slope at TAU(1) is to equal the ! input value C(2,1). ! IBCBEG = 2 means the second derivative at TAU(1) is ! to equal C(2,1). ! IBCEND = 0, 1, or 2 has analogous meaning concerning the ! boundary condition at TAU(N), with the additional ! information taken from C(2,N). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) c(4,n) real ( kind = rk ) divdf1 real ( kind = rk ) divdf3 real ( kind = rk ) dtau real ( kind = rk ) g integer i integer ibcbeg integer ibcend real ( kind = rk ) tau(n) ! ! C(3,*) and C(4,*) are used initially for temporary storage. ! ! Store first differences of the TAU sequence in C(3,*). ! ! Store first divided difference of data in C(4,*). ! do i = 2, n c(3,i) = tau(i) - tau(i-1) end do do i = 2, n c(4,i) = ( c(1,i) - c(1,i-1) ) / ( tau(i) - tau(i-1) ) end do ! ! Construct the first equation from the boundary condition ! at the left endpoint, of the form: ! ! C(4,1) * S(1) + C(3,1) * S(2) = C(2,1) ! ! IBCBEG = 0: Not-a-knot ! if ( ibcbeg == 0 ) then if ( n <= 2 ) then c(4,1) = 1.0D+00 c(3,1) = 1.0D+00 c(2,1) = 2.0D+00 * c(4,2) go to 120 end if c(4,1) = c(3,3) c(3,1) = c(3,2) + c(3,3) c(2,1) = ( ( c(3,2) + 2.0D+00 * c(3,1) ) * c(4,2) * c(3,3) & + c(3,2)**2 * c(4,3) ) / c(3,1) ! ! IBCBEG = 1: derivative specified. ! else if ( ibcbeg == 1 ) then c(4,1) = 1.0D+00 c(3,1) = 0.0D+00 if ( n == 2 ) then go to 120 end if ! ! Second derivative prescribed at left end. ! else c(4,1) = 2.0D+00 c(3,1) = 1.0D+00 c(2,1) = 3.0D+00 * c(4,2) - c(3,2) / 2.0D+00 * c(2,1) if ( n == 2 ) then go to 120 end if end if ! ! If there are interior knots, generate the corresponding ! equations and carry out the forward pass of Gauss elimination, ! after which the I-th equation reads: ! ! C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I). ! do i = 2, n-1 g = -c(3,i+1) / c(4,i-1) c(2,i) = g * c(2,i-1) + 3.0D+00 * ( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) ) c(4,i) = g * c(3,i-1) + 2.0D+00 * ( c(3,i) + c(3,i+1)) end do ! ! Construct the last equation from the second boundary condition, of ! the form ! ! -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N) ! ! If slope is prescribed at right end, one can go directly to ! back-substitution, since the C array happens to be set up just ! right for it at this point. ! if ( ibcend == 1 ) then go to 160 end if if ( 1 < ibcend ) then go to 110 end if ! ! Not-a-knot and 3 <= N, and either 3 < N or also not-a-knot ! at left end point. ! if ( n /= 3 .or. ibcbeg /= 0 ) then g = c(3,n-1) + c(3,n) c(2,n) = ( ( c(3,n) + 2.0D+00 * g ) * c(4,n) * c(3,n-1) + c(3,n)**2 & * ( c(1,n-1) - c(1,n-2) ) / c(3,n-1) ) / g g = - g / c(4,n-1) c(4,n) = c(3,n-1) c(4,n) = c(4,n) + g * c(3,n-1) c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) go to 160 end if ! ! N = 3 and not-a-knot also at left. ! c(2,n) = 2.0D+00 * c(4,n) c(4,n) = 1.0D+00 g = -1.0D+00 / c(4,n-1) c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) go to 160 ! ! IBCEND = 2: Second derivative prescribed at right endpoint. ! 110 continue c(2,n) = 3.0D+00 * c(4,n) + c(3,n) / 2.0D+00 * c(2,n) c(4,n) = 2.0D+00 g = -1.0D+00 / c(4,n-1) c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) go to 160 ! ! N = 2. ! 120 continue if ( ibcend == 2 ) then c(2,n) = 3.0D+00 * c(4,n) + c(3,n) / 2.0D+00 * c(2,n) c(4,n) = 2.0D+00 g = -1.0D+00 / c(4,n-1) c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) else if ( ibcend == 0 .and. ibcbeg /= 0 ) then c(2,n) = 2.0D+00 * c(4,n) c(4,n) = 1.0D+00 g = -1.0D+00 / c(4,n-1) c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) else if ( ibcend == 0 .and. ibcbeg == 0 ) then c(2,n) = c(4,n) end if ! ! Back solve the upper triangular system ! ! C(4,I) * S(I) + C(3,I) * S(I+1) = B(I) ! ! for the slopes C(2,I), given that S(N) is already known. ! 160 continue do i = n-1, 1, -1 c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i) end do ! ! Generate cubic coefficients in each interval, that is, the ! derivatives at its left endpoint, from value and slope at its ! endpoints. ! do i = 2, n dtau = c(3,i) divdf1 = ( c(1,i) - c(1,i-1) ) / dtau divdf3 = c(2,i-1) + c(2,i) - 2.0D+00 * divdf1 c(3,i-1) = 2.0D+00 * ( divdf1 - c(2,i-1) - divdf3 ) / dtau c(4,i-1) = 6.0D+00 * divdf3 / dtau**2 end do return end subroutine cwidth ( w, b, nequ, ncols, integs, nbloks, d, x, iflag ) !*****************************************************************************80 ! !! CWIDTH solves an almost block diagonal linear system. ! ! Discussion: ! ! This routine is a variation of the theme in the algorithm ! by Martin and Wilkinson. It solves the linear system ! A * X = B ! of NEQU equations in case A is almost block diagonal with all ! blocks having NCOLS columns using no more storage than it takes to ! store the interesting part of A. Such systems occur in the determination ! of the B-spline coefficients of a spline approximation. ! ! The block structure of A: ! ! The interesting part of A is taken to consist of NBLOKS ! consecutive blocks, with the I-th block made up of NROWI = INTEGS(1,I) ! consecutive rows and NCOLS consecutive columns of A, and with ! the first LASTI = INTEGS(2,I) columns to the left of the next block. ! These blocks are stored consecutively in the work array W. ! ! For example, here is an 11th order matrix and its arrangement in ! the work array W. (The interesting entries of A are indicated by ! their row and column index modulo 10.) ! ! --- A --- --- W --- ! ! NROW1=3 ! 11 12 13 14 11 12 13 14 ! 21 22 23 24 21 22 23 24 ! 31 32 33 34 NROW2=2 31 32 33 34 ! LAST1=2 43 44 45 46 43 44 45 46 ! 53 54 55 56 NROW3=3 53 54 55 56 ! LAST2=3 66 67 68 69 66 67 68 69 ! 76 77 78 79 76 77 78 79 ! 86 87 88 89 NROW4=1 86 87 88 89 ! LAST3=1 97 98 99 90 NROW5=2 97 98 99 90 ! LAST4=1 08 09 00 01 08 09 00 01 ! 18 19 10 11 18 19 10 11 ! LAST5=4 ! ! For this interpretation of A as an almost block diagonal matrix, ! we have NBLOKS = 5, and the INTEGS array is ! ! I = 1 2 3 4 5 ! K = ! INTEGS(K,I) = 1 3 2 3 1 2 ! 2 2 3 1 1 4 ! ! ! Method: ! ! Gauss elimination with scaled partial pivoting is used, but ! multipliers are not saved in order to save storage. Rather, the ! right hand side is operated on during elimination. The two parameters ! IPVTEQ and LASTEQ are used to keep track of the action. IPVTEQ ! is the index of the variable to be eliminated next, from equations ! IPVTEQ+1,...,LASTEQ, using equation IPVTEQ, possibly after an ! interchange, as the pivot equation. ! ! The entries in the pivot column are always in column ! 1 of W. This is accomplished by putting the entries in rows ! IPVTEQ+1,...,LASTEQ revised by the elimination of the IPVTEQ-th ! variable one to the left in W. In this way, the columns of the ! equations in a given block, as stored in W, will be aligned with ! those of the next block at the moment when these next equations ! become involved in the elimination process. ! ! Thus, for the above example, the first elimination steps proceed ! as follows. ! ! *11 12 13 14 11 12 13 14 11 12 13 14 11 12 13 14 ! *21 22 23 24 *22 23 24 22 23 24 22 23 24 ! *31 32 33 34 *32 33 34 *33 34 33 34 ! 43 44 45 46 43 44 45 46 *43 44 45 46 *44 45 46 ! 53 54 55 56 53 54 55 56 *53 54 55 56 *54 55 56 ! 66 67 68 69 66 67 68 69 66 67 68 69 66 67 68 69 ! ! In all other respects, the procedure is standard, including the ! scaled partial pivoting. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Roger Martin, James Wilkinson, ! Solution of Symmetric and Unsymmetric Band Equations and ! the Calculation of Eigenvectors of Band Matrices, ! Numerische Mathematik, ! Volume 9, Number 4, December 1976, pages 279-301. ! ! Parameters: ! ! Input/output, real ( kind = rk ) W(NEQU,NCOLS), on input, contains ! the interesting part of the almost block diagonal coefficient matrix ! A. The array INTEGS describes the storage scheme. On output, W ! contains the upper triangular factor U of the LU factorization of a ! possibly permuted version of A. In particular, the determinant of ! A could now be found as ! IFLAG * W(1,1) * W(2,1) * ... * W(NEQU,1). ! ! Input/output, real ( kind = rk ) B(NEQU); on input, the right hand ! side of the linear system. On output, B has been overwritten by ! other information. ! ! Input, integer NEQU, the number of equations. ! ! Input, integer NCOLS, the block width, that is, the number of ! columns in each block. ! ! Input, integer INTEGS(2,NEQU), describes the block structure ! of A. ! INTEGS(1,I) = number of rows in block I = NROW. ! INTEGS(2,I) = number of elimination steps in block I = overhang over ! next block = LAST. ! ! Input, integer NBOKS, the number of blocks. ! ! Workspace, real D(NEQU), used to contain row sizes. If storage is ! scarce, the array X could be used in the calling sequence for D. ! ! Output, real ( kind = rk ) X(NEQU), the computed solution, if ! IFLAG is nonzero. ! ! Output, integer IFLAG, error flag. ! = (-1)^(number of interchanges during elimination) if A is invertible; ! = 0 if A is singular. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbloks integer ncols integer nequ real ( kind = rk ) awi1od real ( kind = rk ) b(nequ) real ( kind = rk ) colmax real ( kind = rk ) d(nequ) integer i integer icount integer iflag integer ii integer integs(2,nbloks) integer ipvteq integer ipvtp1 integer istar integer j integer jmax integer lastcl integer lasteq integer lasti integer nexteq integer nrowad real ( kind = rk ) ratio real ( kind = rk ) rowmax real ( kind = rk ) sum1 real ( kind = rk ) temp real ( kind = rk ) w(nequ,ncols) real ( kind = rk ) x(nequ) iflag = 1 ipvteq = 0 lasteq = 0 ! ! The I loop runs over the blocks. ! do i = 1, nbloks ! ! The equations for the current block are added to those currently ! involved in the elimination process, by increasing LASTEQ ! by INTEGS(1,I) after the row size of these equations has been ! recorded in the array D. ! nrowad = integs(1,i) do icount = 1, nrowad nexteq = lasteq + icount rowmax = maxval ( abs ( w(nexteq,1:ncols) ) ) if ( rowmax == 0.0D+00 ) then iflag = 0 return end if d(nexteq) = rowmax end do lasteq = lasteq + nrowad ! ! There will be LASTI = INTEGS(2,I) elimination steps before ! the equations in the next block become involved. ! ! Further, LASTCL records the number of columns involved in the current ! elimination step. It starts equal to NCOLS when a block ! first becomes involved and then drops by one after each elimination ! step. ! lastcl = ncols lasti = integs(2,i) do icount = 1, lasti ipvteq = ipvteq + 1 if ( lasteq <= ipvteq ) then if ( d(ipvteq) < abs ( w(ipvteq,1) ) + d(ipvteq) ) then exit end if iflag = 0 return end if ! ! Determine the smallest ISTAR in (IPVTEQ,LASTEQ) for ! which abs ( W(ISTAR,1) ) / D(ISTAR) is as large as possible, and ! interchange equations IPVTEQ and ISTAR in case IPVTEQ < ISTAR. ! colmax = abs ( w(ipvteq,1) ) / d(ipvteq) istar = ipvteq ipvtp1 = ipvteq + 1 do ii = ipvtp1, lasteq awi1od = abs ( w(ii,1) ) / d(ii) if ( colmax < awi1od ) then colmax = awi1od istar = ii end if end do if ( abs ( w(istar,1) ) + d(istar) == d(istar) ) then iflag = 0 return end if ! ! Rearrange data because of pivoting. ! if ( istar /= ipvteq ) then iflag = -iflag temp = d(istar) d(istar) = d(ipvteq) d(ipvteq) = temp temp = b(istar) b(istar) = b(ipvteq) b(ipvteq) = temp do j = 1, lastcl temp = w(istar,j) w(istar,j) = w(ipvteq,j) w(ipvteq,j) = temp end do end if ! ! Subtract the appropriate multiple of equation IPVTEQ from ! equations IPVTEQ+1,...,LASTEQ to make the coefficient of the ! IPVTEQ-th unknown (presently in column 1 of W) zero, but ! store the new coefficients in W one to the left from the old. ! do ii = ipvtp1, lasteq ratio = w(ii,1) / w(ipvteq,1) do j = 2, lastcl w(ii,j-1) = w(ii,j) - ratio * w(ipvteq,j) end do w(ii,lastcl) = 0.0D+00 b(ii) = b(ii) - ratio * b(ipvteq) end do lastcl = lastcl - 1 end do end do ! ! At this point, W and B contain an upper triangular linear system ! equivalent to the original one, with W(I,J) containing entry ! (I, I-1+J) of the coefficient matrix. Solve this system by ! back substitution, taking into account its block structure. ! ! I-loop over the blocks, in reverse order. ! i = nbloks do while ( 0 < i ) lasti = integs(2,i) jmax = ncols - lasti do icount = 1, lasti sum1 = dot_product ( x(ipvteq+1:ipvteq+jmax), w(ipvteq,2:jmax+1) ) x(ipvteq) = ( b(ipvteq) - sum1 ) / w(ipvteq,1) jmax = jmax + 1 ipvteq = ipvteq - 1 end do i = i - 1 end do return end subroutine difequ ( mode, xx, v ) !*****************************************************************************80 ! !! DIFEQU returns information about a differential equation. ! ! Discussion: ! ! This sample version of DIFEQU is for the example in chapter XV. It is a ! nonlinear second order two point boundary value problem. ! ! Modified: ! ! 16 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, integer MODE, an integer indicating the task to ! be performed. ! 1, initialization ! 2, evaluate the differential equation at point XX. ! 3, specify the next side condition ! 4, analyze the approximation ! ! Input, real ( kind = rk ) XX, a point at which information is wanted ! ! Output, real ( kind = rk ) V, depends on the MODE. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: npiece = 100 integer, parameter :: ncoef = 2000 real ( kind = rk ) break real ( kind = rk ) coef real ( kind = rk ), save :: eps real ( kind = rk ) ep1 real ( kind = rk ) ep2 real ( kind = rk ) error real ( kind = rk ), save :: factor integer i integer iside integer itermx integer k integer kpm integer l integer m integer mode real ( kind = rk ) ppvalu real ( kind = rk ) rho real ( kind = rk ), save :: s2ovep real ( kind = rk ) solutn real ( kind = rk ) un real ( kind = rk ) v(20) real ( kind = rk ) value real ( kind = rk ) x real ( kind = rk ) xside real ( kind = rk ) xx save / approx / save / other / save / side / common / approx / break(npiece), coef(ncoef), l, kpm common / other / itermx, k, rho(19) common / side / m, iside, xside(10) ! ! Initialize everything, Set the order M of the differential equation, ! the nondecreasing sequence XSIDE(1:M), of points at which side ! conditions are given and anything else necessary. ! if ( mode == 1 ) then m = 2 xside(1) = 0.0D+00 xside(2) = 1.0D+00 ! ! Print out heading. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Carrier''s nonlinear perturbation problem' write ( *, '(a)' ) ' ' eps = 0.005D+00 write ( *, '(a,g14.6)' ) ' EPS = ', eps ! ! Set constants used in formula for solution below. ! factor = ( sqrt ( 2.0D+00 ) + sqrt ( 3.0D+00 ) )**2 s2ovep = sqrt ( 2.0D+00 / eps ) ! ! Initial guess for Newton iteration: UN(X) = X*X-1. ! l = 1 break(1) = 0.0D+00 coef(1:kpm) = 0.0D+00 coef(1) = -1.0D+00 coef(3) = 2.0D+00 itermx = 10 ! ! Provide value of left side coefficients and right hand side at XX. ! Specifically, at XX the differential equation reads: ! ! V(M+1) D^M + V(M) D^(M-1) + ... + V(1) D^0 = V(M+2) ! ! in terms of the quantities V(1:M+2), to be computed here. ! else if ( mode == 2 ) then v(3) = eps v(2) = 0.0D+00 un = ppvalu ( break, coef, l, kpm, xx, 0 ) v(1) = 2.0D+00 * un v(4) = un**2 + 1.0D+00 ! ! Provide the M side conditions. these conditions are of the form ! ! V(M+1) D^M + V(M) D^(M-1) + ... + V(1) D^0 = V(M+2) ! ! in terms of the quantities V(1:M+2), to be specified here. ! Note that V(M+1) = 0 for customary side conditions. ! else if ( mode == 3 ) then v(m+1) = 0.0D+00 if ( iside == 1 ) then v(2) = 1.0D+00 v(1) = 0.0D+00 v(4) = 0.0D+00 iside = iside + 1 else if ( iside == 2 ) then v(2) = 0.0D+00 v(1) = 1.0D+00 v(4) = 0.0D+00 iside = iside + 1 end if ! ! Calculate the error near the boundary layer at 1. ! else if ( mode == 4 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' X G(X) G(X)-F(X):' write ( *, '(a)' ) ' ' x = 0.75D+00 do i = 1, 9 ep1 = exp ( s2ovep * ( 1.0D+00 - x ) ) * factor ep2 = exp ( s2ovep * ( 1.0D+00 + x ) ) * factor solutn = 12.0D+00 / ( 1.0D+00 + ep1 )**2 * ep1 & + 12.0D+00 / ( 1.0D+00 + ep2 )**2 * ep2 - 1.0D+00 value = ppvalu ( break, coef, l, kpm, x, 0 ) error = solutn - value write ( *, '(2x,3g14.6)' ) x, solutn, error x = x + 0.03125D+00 end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIFEQU - Fatal error!' write ( *, '(a)' ) ' Illegal value of MODE:' write ( *, '(a,i8)' ) mode stop 1 end if return end subroutine dtblok ( bloks, integs, nbloks, ipivot, iflag, detsgn, detlog ) !*****************************************************************************80 ! !! DTBLOK gets the determinant of an almost block diagonal matrix. ! ! Discussion: ! ! The matrix's PLU factorization must have been obtained ! previously by FCBLOK. ! ! The logarithm of the determinant is computed instead of the ! determinant itself to avoid the danger of overflow or underflow ! inherent in this calculation. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BLOKS(*), the factorization of A computed ! by FCBLOK. ! ! Input, integer INTEGS(3,NBLOKS), describes the block ! structure of A. ! ! Input, integer NBLOKS, the number of blocks in A. ! ! Input, integer IPIVOT(*), pivoting information. ! The dimension of IPIVOT is the sum ( INTEGS(1,1:NBLOKS) ). ! ! Input, integer IFLAG, = (-1)^(number of interchanges during ! factorization) if successful, otherwise IFLAG = 0. ! ! Output, real ( kind = rk ) DETSGN, the sign of the determinant. ! ! Output, real ( kind = rk ) DETLOG, the natural logarithm of the ! determinant, if the determinant is not zero. If the determinant ! is 0, then DETLOG is returned as 0. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbloks real ( kind = rk ) bloks(*) real ( kind = rk ) detlog real ( kind = rk ) detsgn integer i integer iflag integer index integer indexp integer integs(3,nbloks) integer ip integer ipivot(1) integer k integer last integer nrow detsgn = iflag detlog = 0.0D+00 if ( iflag == 0 ) then return end if index = 0 indexp = 0 do i = 1, nbloks nrow = integs(1,i) last = integs(3,i) do k = 1, last ip = index + nrow * ( k - 1 ) + ipivot(indexp+k) detlog = detlog + log ( abs ( bloks(ip) ) ) detsgn = detsgn * sign ( 1.0D+00, bloks(ip) ) end do index = nrow * integs(2,i) + index indexp = indexp + nrow end do return end subroutine eqblok ( t, n, kpm, work1, work2, bloks, lenblk, integs, nbloks, b ) !*****************************************************************************80 ! !! EQBLOK is to be called in COLLOC. ! ! Method: ! ! Each breakpoint interval gives rise to a block in the linear system. ! This block is determined by the K collocation equations in the interval ! with the side conditions, if any, in the interval interspersed ! appropriately, and involves the KPM B-splines having the interval in ! their support. Correspondingly, such a block has NROW = K + ISIDEL ! rows, with ISIDEL = number of side conditions in this and the ! previous intervals, and NCOL = KPM columns. ! ! Further, because the interior knots have multiplicity K, we can ! carry out in SLVBLK K elimination steps in a block before pivoting ! might involve an equation from the next block. In the last block, ! of course, all KPM elimination steps will be carried out in SLVBLK. ! ! See the detailed comments in SLVBLK for further ! information about the almost block diagonal form used here. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(N+KPM), the knot sequence. ! ! Input, integer N, the dimension of the approximating spline ! space, that is, the order of the linear system to be constructed. ! ! Input, integer KPM, = K + M, the order of the approximating ! spline. ! ! Input, integer LENBLK, the maximum length of the array BLOKS, ! as allowed by the dimension statement in COLLOC. ! ! Workspace, real ( kind = rk ) WORK1(KPM,KPM), used in PUTIT. ! ! Workspace, real ( kind = rk ) WORK2(KPM,M+1), used in PUTIT. ! ! Output, real ( kind = rk ) BLOKS(*), the coefficient matrix of the ! linear system, stored in almost block diagonal form, of size ! KPM * sum ( INTEGS(1,1:NBLOKS) ). ! ! Output, integer INTEGS(3,NBLOKS), describing the block ! structure. ! INTEGS(1,I) = number of rows in block I; ! INTEGS(2,I) = number of columns in block I; ! INTEGS(3,I) = number of elimination steps which can be carried out in ! block I before pivoting might bring in an equation from the next block. ! ! Output, integer NBLOKS, the number of blocks, equals number ! of polynomial pieces. ! ! Output, real ( kind = rk ) B(*), the right hand side of the linear ! system, stored corresponding to the almost block diagonal form, ! of size sum ( INTEGS(1,1:NBLOKS) ). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer kpm integer n real ( kind = rk ) b(*) real ( kind = rk ) bloks(*) integer i integer index integer indexb integer integs(3,*) integer iside integer isidel integer itermx integer k integer left integer lenblk integer m integer nbloks integer nrow real ( kind = rk ) rho real ( kind = rk ) t(n+kpm) real ( kind = rk ) work1(kpm,kpm) real ( kind = rk ) work2(kpm,*) real ( kind = rk ) xside save / other / save / side / common / other / itermx, k, rho(19) common / side / m, iside, xside(10) index = 1 indexb = 1 i = 0 iside = 1 do left = kpm, n, k i = i + 1 ! ! Determine INTEGS(:,I). ! integs(2,i) = kpm if ( n <= left ) then integs(3,i) = kpm isidel = m ! ! At this point, ISIDE - 1 gives the number of side conditions ! incorporated so far. Adding to this the side conditions in the ! current interval gives the number ISIDEL. ! else integs(3,i) = k isidel = iside - 1 do if ( isidel == m ) then exit end if if ( t(left+1) <= xside(isidel+1) ) then exit end if isidel = isidel + 1 end do end if nrow = k + isidel integs(1,i) = nrow ! ! The detailed equations for this block are generated and put ! together in PUTIT. ! if ( lenblk < index + nrow * kpm - 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EQBLOK - Fatal error!' write ( *, '(a)' ) ' The dimension of BLOKS is too small.' write ( *, '(a,i8)' ) ' LENBLK = ', lenblk stop 1 end if call putit ( t, kpm, left, work1, work2, bloks(index), nrow, b(indexb) ) index = index + nrow * kpm indexb = indexb + nrow end do nbloks = i return end subroutine evnnot ( break, coef, l, k, brknew, lnew, coefg ) !*****************************************************************************80 ! !! EVNNOT is a version of NEWNOT returning uniform knots. ! ! Discussion: ! ! EVNNOT returns LNEW+1 knots in BRKNEW which are evenly spaced between ! BREAK(1) and BREAK(L+1). ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BREAK(L+1), real ( kind = rk ) COEF(K,L), ! integer L, integer K, the piecewise polynomial representation ! of a certain function F of order K. Specifically, ! d^(K-1) F(X) = COEF(K,I) for BREAK(I) <= X < BREAK(I+1). ! ! Input, integer LNEW, the number of subintervals into which ! the interval (A,B) is to be sectioned by the new breakpoint ! sequence BRKNEW. ! ! Output, real ( kind = rk ) BRKNEW(LNEW+1), the new breakpoints. ! ! Output, real (kind = rk ) COEFG(2,L), the coefficient part of the ! piecewise polynomial representation BREAK, COEFG, L, 2 for the monotone ! piecewise linear function G with respect to which BRKNEW will ! be equidistributed. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer l integer lnew real ( kind = rk ) break(l+1) real ( kind = rk ) brknew(lnew+1) real ( kind = rk ) coef(k,l) real ( kind = rk ) coefg(2,l) integer i coefg(2,l) = 0.0D+00 if ( lnew == 0 ) then brknew(1) = 0.5D+00 * ( break(1) + break(l+1) ) else do i = 1, lnew + 1 brknew(i) = ( real ( lnew - i + 1, kind = rk ) * break(1) & + real ( i - 1, kind = rk ) * break(l+1) ) & / real ( lnew, kind = rk ) end do end if return end subroutine factrb ( w, ipivot, d, nrow, ncol, last, iflag ) !*****************************************************************************80 ! !! FACTRB constructs a partial PLU factorization. ! ! Discussion: ! ! This factorization corresponds to steps 1 through LAST in Gauss ! elimination for the matrix W of order ( NROW, NCOL ), using ! pivoting of scaled rows. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input/output, real ( kind = rk ) W(NROW,NCOL); on input, contains the ! matrix to be partially factored; on output, the partial factorization. ! ! Output, integer IPIVOT(NROW), contains a record of the ! pivoting strategy used; row IPIVOT(I) is used during the I-th elimination ! step, for I = 1, ..., LAST. ! ! Workspace, real ( kind = rk ) D(NROW), used to store the maximum entry ! in each row. ! ! Input, integer NROW, the number of rows of W. ! ! Input, integer NCOL, the number of columns of W. ! ! Input, integer LAST, the number of elimination steps to ! be carried out. ! ! Input/output, integer IFLAG. On output, equals the input ! value times (-1)^(number of row interchanges during the factorization ! process), in case no zero pivot was encountered. ! Otherwise, IFLAG = 0 on output. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ncol integer nrow real ( kind = rk ) awikdi real ( kind = rk ) colmax real ( kind = rk ) d(nrow) integer i integer iflag integer ipivi integer ipivk integer ipivot(nrow) integer j integer k integer kp1 integer last real ( kind = rk ) ratio real ( kind = rk ) rowmax real ( kind = rk ) w(nrow,ncol) ! ! Initialize IPIVOT and D. ! do i = 1, nrow ipivot(i) = i end do do i = 1, nrow rowmax = maxval ( abs ( w(i,1:ncol) ) ) if ( rowmax == 0.0D+00 ) then iflag = 0 return end if d(i) = rowmax end do ! ! Gauss elimination with pivoting of scaled rows, loop over K = 1,..., LAST. ! k = 1 ! ! As pivot row for K-th step, pick among the rows not yet used, ! that is, from rows IPIVOT(K:NROW), the one whose K-th entry, compared ! to the row size, is largest. ! ! If this row does not turn out to be row IPIVOT(K), redefine IPIVOT(K) ! appropriately and record this interchange by changing the sign ! of IFLAG. ! do while ( k <= last ) ipivk = ipivot(k) if ( k == nrow ) then if ( abs ( w(ipivk,nrow) ) + d(ipivk) <= d(ipivk) ) then iflag = 0 end if return end if j = k kp1 = k + 1 colmax = abs ( w(ipivk,k) ) / d(ipivk) ! ! Find the largest pivot. ! do i = kp1, nrow ipivi = ipivot(i) awikdi = abs ( w(ipivi,k) ) / d(ipivi) if ( colmax < awikdi ) then colmax = awikdi j = i end if end do if ( j /= k ) then ipivk = ipivot(j) ipivot(j) = ipivot(k) ipivot(k) = ipivk iflag = - iflag end if ! ! If the pivot element is too small in absolute value, declare ! the matrix to be noninvertible and quit. ! if ( abs ( w(ipivk,k) ) + d(ipivk) <= d(ipivk) ) then iflag = 0 return end if ! ! Otherwise, subtract the appropriate multiple of the pivot ! row from the remaining rows, that is, the rows IPIVOT(K+1:NROW), ! to make the K-th entry zero. ! ! Save the multiplier in its place. ! do i = kp1, nrow ipivi = ipivot(i) w(ipivi,k) = w(ipivi,k) / w(ipivk,k) ratio = - w(ipivi,k) w(ipivi,kp1:ncol) = ratio * w(ipivk,kp1:ncol) + w(ipivi,kp1:ncol) end do k = kp1 end do return end subroutine fcblok ( bloks, integs, nbloks, ipivot, scrtch, iflag ) !*****************************************************************************80 ! !! FCBLOK supervises the PLU factorization of an almost block diagonal matrix. ! ! Discussion: ! ! The routine supervises the PLU factorization with pivoting of ! the scaled rows of an almost block diagonal matrix. ! ! The almost block diagonal matrix is stored in the arrays ! BLOKS and INTEGS. ! ! The FACTRB routine carries out steps 1,..., LAST of Gauss ! elimination, with pivoting, for an individual block. ! ! The SHIFTB routine shifts the remaining rows to the top of ! the next block. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input/output, real ( kind = rk ) BLOKS(*). On input, the almost ! block diagonal matrix A to be factored. On output, the ! factorization of A. ! ! Input, integer INTEGS(3,NBLOKS), describes the block ! structure of A. ! ! Input, integer NBLOKS, the number of blocks in A. ! ! Output, integer IPIVOT(*), which will contain pivoting ! information. The dimension of IPIVOT is the sum ( INTEGS(1,1:NBLOKS) ). ! ! Workspace, real SCRTCH(*), of length maxval ( integs(1,1:NBLOKS) ). ! ! Output, integer IFLAG, error flag. ! = 0, in case matrix was found to be singular; ! = (-1)^(number of row interchanges during factorization), otherwise. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbloks real ( kind = rk ) bloks(*) integer i integer iflag integer index integer indexb integer indexn integer integs(3,nbloks) integer ipivot(*) integer last integer ncol integer nrow real ( kind = rk ) scrtch(*) iflag = 1 indexb = 1 indexn = 1 i = 1 ! ! Loop over the blocks. I is the loop index. ! do index = indexn nrow = integs(1,i) ncol = integs(2,i) last = integs(3,i) ! ! Carry out elimination on the I-th block until next block ! enters, for columns 1 through LAST of I-th block. ! call factrb ( bloks(index), ipivot(indexb), scrtch, nrow, ncol, & last, iflag ) ! ! Check for having reached a singular block or the last block. ! if ( iflag == 0 .or. i == nbloks ) then exit end if i = i + 1 indexn = nrow * ncol + index ! ! Put the rest of the I-th block onto the next block. ! call shiftb ( bloks(index), ipivot(indexb), nrow, ncol, last, & bloks(indexn), integs(1,i), integs(2,i) ) indexb = indexb + nrow end do return end subroutine interv ( xt, lxt, x, left, mflag ) !*****************************************************************************80 ! !! INTERV brackets a real value in an ascending vector of values. ! ! Discussion: ! ! The XT array is a set of increasing values. The goal of the routine ! is to determine the largest index I so that ! ! XT(I) < XT(LXT) and XT(I) <= X. ! ! The routine is designed to be efficient in the common situation ! that it is called repeatedly, with X taken from an increasing ! or decreasing sequence. ! ! This will happen when a piecewise polynomial is to be graphed. ! The first guess for LEFT is therefore taken to be the value ! returned at the previous call and stored in the local variable ILO. ! ! A first check ascertains that ILO < LXT. This is necessary ! since the present call may have nothing to do with the previous ! call. Then, if ! XT(ILO) <= X < XT(ILO+1), ! we set LEFT = ILO and are done after just three comparisons. ! ! Otherwise, we repeatedly double the difference ISTEP = IHI - ILO ! while also moving ILO and IHI in the direction of X, until ! XT(ILO) <= X < XT(IHI) ! after which we use bisection to get, in addition, ILO + 1 = IHI. ! The value LEFT = ILO is then returned. ! ! Thanks to Daniel Gloger for pointing out an important modification ! to the routine, so that the piecewise polynomial in B-form is ! left-continuous at the right endpoint of the basic interval, ! 17 April 2014. ! ! Modified: ! ! 17 April 2014 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) XT(LXT), a nondecreasing sequence of values. ! ! Input, integer LXT, the dimension of XT. ! ! Input, real ( kind = rk ) X, the point whose location with ! respect to the sequence XT is to be determined. ! ! Output, integer LEFT, the index of the bracketing value: ! 1 if X < XT(1) ! I if XT(I) <= X < XT(I+1) ! I if XT(I) < X == XT(I+1) == XT(LXT) ! ! Output, integer MFLAG, indicates whether X lies within the ! range of the data. ! -1: X < XT(1) ! 0: XT(I) <= X < XT(I+1) ! +1: XT(LXT) < X ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer lxt integer left integer mflag integer ihi integer, save :: ilo = 1 integer istep integer middle real ( kind = rk ) x real ( kind = rk ) xt(lxt) ihi = ilo + 1 if ( lxt <= ihi ) then if ( xt(lxt) <= x ) then go to 110 end if if ( lxt <= 1 ) then mflag = -1 left = 1 return end if ilo = lxt - 1 ihi = lxt end if if ( xt(ihi) <= x ) then go to 20 end if if ( xt(ilo) <= x ) then mflag = 0 left = ilo return end if ! ! Now X < XT(ILO). Decrease ILO to capture X. ! istep = 1 10 continue ihi = ilo ilo = ihi - istep if ( 1 < ilo ) then if ( xt(ilo) <= x ) then go to 50 end if istep = istep * 2 go to 10 end if ilo = 1 if ( x < xt(1) ) then mflag = -1 left = 1 return end if go to 50 ! ! Now XT(IHI) <= X. Increase IHI to capture X. ! 20 continue istep = 1 30 continue ilo = ihi ihi = ilo + istep if ( ihi < lxt ) then if ( x < xt(ihi) ) then go to 50 end if istep = istep * 2 go to 30 end if if ( xt(lxt) <= x ) then go to 110 end if ! ! Now XT(ILO) < = X < XT(IHI). Narrow the interval. ! ihi = lxt 50 continue do middle = ( ilo + ihi ) / 2 if ( middle == ilo ) then mflag = 0 left = ilo return end if ! ! It is assumed that MIDDLE = ILO in case IHI = ILO+1. ! if ( xt(middle) <= x ) then ilo = middle else ihi = middle end if end do ! ! Set output and return. ! 110 continue mflag = 1 if ( x == xt(lxt) ) then mflag = 0 end if do left = lxt - 1, 1, -1 if ( xt(left) < xt(lxt) ) then return end if end do return end subroutine knots ( break, l, kpm, m, t, n ) !*****************************************************************************80 ! !! KNOTS is to be called in COLLOC. ! ! Discussion: ! ! Note that the FORTRAN77 calling sequence has been modified, by ! adding the variable M. ! ! From the given breakpoint sequence BREAK, this routine constructs the ! knot sequence T so that ! SPLINE(K+M,T) = PP(K+M,BREAK) ! with M-1 continuous derivatives. ! ! This means that T(1:N+KPM) is equal to BREAK(1) KPM times, then ! BREAK(2) through BREAK(L) each K times, then, finally, BREAK(L+1) ! KPM times. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BREAK(L+1), the breakpoint sequence. ! ! Input, integer L, the number of intervals or pieces. ! ! Input, integer KPM, = K+M, the order of the piecewise ! polynomial function or spline. ! ! Input, integer M, the order of the differential equation. ! ! Output, real ( kind = rk ) T(N+KPM), the knot sequence. ! ! Output, integer N, = L*K+M = the dimension of SPLINE(K+M,T). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer kpm integer l integer n real ( kind = rk ) break(l+1) integer j integer jj integer jjj integer k integer ll integer m real ( kind = rk ) t(n+kpm) k = kpm - m n = l * k + m jj = n + kpm jjj = l + 1 do ll = 1, kpm t(jj) = break(jjj) jj = jj - 1 end do do j = 1, l jjj = jjj - 1 do ll = 1, k t(jj) = break(jjj) jj = jj - 1 end do end do t(1:kpm) = break(1) return end subroutine l2appr ( t, n, k, q, diag, bcoef ) !*****************************************************************************80 ! !! L2APPR constructs a weighted L2 spline approximation to given data. ! ! Discussion: ! ! The routine constructs the weighted discrete L2-approximation by ! splines of order K with knot sequence T(1:N+K) to ! given data points ( TAU(1:NTAU), GTAU(1:NTAU) ). ! ! The B-spline coefficients BCOEF of the approximating spline are ! determined from the normal equations using Cholesky's method. ! ! Method: ! ! The B-spline coefficients of the L2-approximation are determined as the ! solution of the normal equations, for 1 <= I <= N: ! sum ( 1 <= J <= N ) ( B(I), B(J) ) * BCOEF(J) = ( B(I), G ). ! ! Here, B(I) denotes the I-th B-spline, G denotes the function to ! be approximated, and the inner product of two functions F and G ! is given by ! ( F, G ) = sum ( 1 <= I <= NTAU ) WEIGHT(I) * F(TAU(I)) * G(TAU(I)). ! ! The arrays TAU and WEIGHT are given in common block DATA, as is the ! array GTAU(1:NTAU) = G(TAU(1:NTAU)). ! ! The values of the B-splines B(1:N) are supplied by BSPLVB. ! ! The coefficient matrix C, with ! C(I,J) = ( B(I), B(J) ) ! of the normal equations is symmetric and (2*K-1)-banded, therefore ! can be specified by giving its K bands at or below the diagonal. ! ! For I = 1:N and J = I:min(I+K-1,N), we store ! ( B(I), B(J) ) = C(I,J) ! in ! Q(I-J+1,J), ! and the right hand side ! ( B(I), G ) ! in ! BCOEF(I). ! ! Since B-spline values are most efficiently generated by finding ! simultaneously the value of every nonzero B-spline at one point, ! the entries of C (that is, of Q), are generated by computing, for ! each LL, all the terms involving TAU(LL) simultaneously and adding ! them to all relevant entries. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. ! ! Input, integer N, the dimension of the space of splines ! of order K with knots T. ! ! Input, integer K, the order of the splines. ! ! Workspace, real ( kind = rk ) Q(K,N), used to store the K lower ! diagonals of the Gramian matrix C. ! ! Workspace, real ( kind = rk ) DIAG(N), used in BCHFAC. ! ! Output, real ( kind = rk ) BCOEF(N), the B-spline coefficients of ! the L2 approximation to the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer n integer, parameter :: ntmax = 200 real ( kind = rk ) bcoef(n) real ( kind = rk ) biatx(k) real ( kind = rk ) diag(n) real ( kind = rk ) dw real ( kind = rk ) gtau integer i integer j integer jj integer left integer leftmk integer ll integer mm integer ntau real ( kind = rk ) q(k,n) real ( kind = rk ) t(n+k) real ( kind = rk ) tau real ( kind = rk ) totalw real ( kind = rk ) weight save / i4data / save / r8data / common / i4data / ntau common / r8data / tau(ntmax), gtau(ntmax), weight(ntmax), totalw bcoef(1:n) = 0.0D+00 q(1:k,1:n) = 0.0D+00 left = k leftmk = 0 do ll = 1, ntau ! ! Locate LEFT such that TAU(LL) is in ( T(LEFT), T(LEFT+1) ). ! do if ( left == n ) then exit end if if ( tau(ll) < t(left+1) ) then exit end if left = left + 1 leftmk = leftmk + 1 end do call bsplvb ( t, k, 1, tau(ll), left, biatx ) ! ! BIATX(MM) contains the value of B(LEFT-K+MM) at TAU(LL). ! ! Hence, with DW = BIATX(MM) * WEIGHT(LL), the number DW * GTAU(LL) ! is a summand in the inner product ! ( B(LEFT-K+MM), G) ! which goes into BCOEF(LEFT-K+MM) ! and the number BIATX(JJ)*DW is a summand in the inner product ! (B(LEFT-K+JJ), B(LEFT-K+MM)), into Q(JJ-MM+1,LEFT-K+MM) ! since (LEFT-K+JJ)-(LEFT-K+MM)+1 = JJ - MM + 1. ! do mm = 1, k dw = biatx(mm) * weight(ll) j = leftmk + mm bcoef(j) = dw * gtau(ll) + bcoef(j) i = 1 do jj = mm, k q(i,j) = biatx(jj) * dw + q(i,j) i = i + 1 end do end do end do ! ! Construct the Cholesky factorization for C in Q, then ! use it to solve the normal equations ! ! C * X = BCOEF ! ! for X, and store X in BCOEF. ! call bchfac ( q, k, n, diag ) call bchslv ( q, k, n, bcoef ) return end subroutine l2err ( iprfun, ftau, error ) !*****************************************************************************80 ! !! L2ERR computes the errors of an L2 approximation. ! ! Discussion: ! ! This routine computes various errors of the current L2 approximation, ! whose piecewise polynomial representation is contained in common ! block APPROX, to the given data contained in common block DATA. ! ! It prints out the average error ERRL1, the L2 error ERRL2, and the ! maximum error ERRMAX. ! ! Modified: ! ! 16 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, integer IPRFUN. If IPRFUN = 1, the routine prints out ! the value of the approximation as well as its error at ! every data point. ! ! Output, real ( kind = rk ) FTAU(NTAU), contains the value of the computed ! approximation at each value TAU(1:NTAU). ! ! Output, real ( kind = rk ) ERROR(NTAU), with ! ERROR(I) = SCALE * ( G - F )(TAU(I)). Here, SCALE equals 1 ! in case IPRFUN /= 1, or the absolute error is greater than 100 ! somewhere. Otherwise, SCALE is such that the maximum of the ! absolute value of ERROR(1:NTAU) lies between 10 and 100. This ! makes the printed output more illustrative. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: lpkmax = 100 integer, parameter :: ntmax = 200 integer, parameter :: ltkmax = 2000 integer ntau real ( kind = rk ) break real ( kind = rk ) coef real ( kind = rk ) err real ( kind = rk ) errl1 real ( kind = rk ) errl2 real ( kind = rk ) errmax real ( kind = rk ) error(ntau) real ( kind = rk ) ftau(ntau) real ( kind = rk ) gtau integer ie integer iprfun integer k integer l integer ll real ( kind = rk ) ppvalu real ( kind = rk ) scale real ( kind = rk ) tau real ( kind = rk ) totalw real ( kind = rk ) weight save / approx / save / i4data / save / r8data / common / approx / break(lpkmax), coef(ltkmax), l, k common / i4data / ntau common / r8data / tau(ntmax), gtau(ntmax), weight(ntmax), totalw errl1 = 0.0D+00 errl2 = 0.0D+00 errmax = 0.0D+00 do ll = 1, ntau ftau(ll) = ppvalu ( break, coef, l, k, tau(ll), 0 ) error(ll) = gtau(ll) - ftau(ll) err = abs(error(ll)) if ( errmax < err ) then errmax = err end if errl1 = errl1 + err * weight(ll) errl2 = errl2 + err**2 * weight(ll) end do errl1 = errl1 / totalw errl2 = sqrt ( errl2 / totalw ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Least square error = ', errl2 write ( *, '(a,g14.6)' ) ' Average error = ', errl1 write ( *, '(a,g14.6)' ) ' Maximum error = ', errmax write ( *, '(a)' ) ' ' if ( iprfun /= 1 ) then return end if ! ! Scale error curve and print. ! ie = 0 scale = 1.0D+00 if ( errmax < 10.0D+00 ) then do ie = 1, 9 scale = scale * 10.0D+00 if ( 10.0D+00 <= errmax * scale ) then exit end if end do end if error(1:ntau) = error(1:ntau) * scale write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Approximation and scaled error curve' write ( *, '(a)' ) ' ' write ( *, '(a,i1)' ) & ' Data point Approximation Deviation x 10**', ie write ( *, '(a)' ) ' ' write ( *, '(i4,f16.8,f16.8,f17.6)' ) & ( ll, tau(ll), ftau(ll), error(ll), ll = 1, ntau ) return end subroutine l2knts ( break, l, k, t, n ) !*****************************************************************************80 ! !! L2KNTS converts breakpoints to knots. ! ! Discussion: ! ! The breakpoint sequence BREAK is converted into a corresponding ! knot sequence T to allow the representation of a piecewise ! polynomial function of order K with K-2 continuous derivatives ! as a spline of order K with knot sequence T. ! ! This means that T(1:N+K) = BREAK(1) K times, then BREAK(2:L), ! then BREAK(L+1) K times. ! ! Therefore, N = K - 1 + L. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, integer K, the order. ! ! Input, integer L, the number of polynomial pieces. ! ! Input, real ( kind = rk ) BREAK(L+1), the breakpoint sequence. ! ! Output, real ( kind = rk ) T(N+K), the knot sequence. ! ! Output, integer N, the dimension of the corresponding spline ! space of order K. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer l integer n real ( kind = rk ) break(l+1) real ( kind = rk ) t(k-1+l+k) n = k - 1 + l t(1:k-1) = break(1) t(k:n) = break(1:l) t(n+1:n+k) = break(l+1) return end subroutine newnot ( break, coef, l, k, brknew, lnew, coefg ) !*****************************************************************************80 ! !! NEWNOT returns LNEW+1 knots which are equidistributed on (A,B). ! ! Discussion: ! ! The knots are equidistributed on (A,B) = ( BREAK(1), BREAK(L+1) ) ! with respect to a certain monotone function G related to the K-th root of ! the K-th derivative of the piecewise polynomial function F whose ! piecewise polynomial representation is contained in BREAK, COEF, L, K. ! ! Method: ! ! The K-th derivative of the given piecewise polynomial function F does ! not exist, except perhaps as a linear combination of delta functions. ! ! Nevertheless, we construct a piecewise constant function H with ! breakpoint sequence BREAK which is approximately proportional ! to abs ( d^K(F) ). ! ! Specifically, on (BREAK(I), BREAK(I+1)), ! ! abs(jump at BREAK(I) of PC) abs(jump at BREAK(I+1) of PC) ! H = --------------------------- + ---------------------------- ! BREAK(I+1) - BREAK(I-1) BREAK(I+2) - BREAK(I) ! ! with PC the piecewise constant (K-1)st derivative of F. ! ! Then, the piecewise linear function G is constructed as ! G(X) = integral ( A <= Y <= X ) H(Y)^(1/K) dY, ! and its piecewise polynomial coefficients are stored in COEFG. ! ! Then BRKNEW is determined by ! BRKNEW(I) = A + G^(-1)((I-1)*STEP), for I = 1:LNEW+1, ! where STEP = G(B) / LNEW and (A,B) = ( BREAK(1), BREAK(L+1) ). ! ! In the event that PC = d^(K-1)(F) is constant in ( A, B ) and ! therefore H = 0 identically, BRKNEW is chosen uniformly spaced. ! ! If IPRINT is set positive, then the piecewise polynomial coefficients ! of G will be printed out. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BREAK(L+1), real ( kind = rk ) COEF(K,L), ! integer L, integer K, the piecewise polynomial representation ! of a certain function F of order K. Specifically, ! d^(k-1) F(X) = COEF(K,I) for BREAK(I) <= X < BREAK(I+1). ! ! Input, integer LNEW, the number of intervals into which the ! interval (A,B) is to be divided by the new breakpoint sequence BRKNEW. ! ! Output, real ( kind = rk ) BRKNEW(LNEW+1), the new breakpoint sequence. ! ! Output, real ( kind = rk ) COEFG(2,L), the coefficient part of the piecewise ! polynomial representation BREAK, COEFG, L, 2 for the monotone piecewise ! linear function G with respect to which BRKNEW will be equidistributed. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer l integer lnew real ( kind = rk ) break(l+1) real ( kind = rk ) brknew(lnew+1) real ( kind = rk ) coef(k,l) real ( kind = rk ) coefg(2,l) real ( kind = rk ) dif real ( kind = rk ) difprv integer i integer, save :: iprint = 0 integer j real ( kind = rk ) oneovk real ( kind = rk ) step real ( kind = rk ) stepi ! ! If G is constant, BRKNEW is uniform. ! if ( l <= 1 ) then call evnnot ( break, coef, l, k, brknew, lnew, coefg ) return end if brknew(1) = break(1) brknew(lnew+1) = break(l+1) ! ! Construct the continuous piecewise linear function G. ! oneovk = 1.0D+00 / real ( k, kind = rk ) coefg(1,1) = 0.0D+00 difprv = abs ( coef(k,2) - coef(k,1) ) / ( break(3) - break(1) ) do i = 2, l dif = abs ( coef(k,i) - coef(k,i-1) ) / ( break(i+1) - break(i-1) ) coefg(2,i-1) = ( dif + difprv )**oneovk coefg(1,i) = coefg(1,i-1) + coefg(2,i-1) * ( break(i) - break(i-1) ) difprv = dif end do coefg(2,l) = ( 2.0D+00 * difprv )**oneovk ! ! STEP = G(B) / LNEW. ! step = ( coefg(1,l) + coefg(2,l) * ( break(l+1) - break(l) ) ) & / real ( lnew, kind = rk ) if ( 0 < iprint ) then write ( *, '(2x,e16.7)' ) step do i = 1, l write ( *, '(i5,2e16.5)' ) i, coefg(1:2,i) end do end if ! ! If G is constant, BRKNEW is uniform. ! if ( step <= 0.0D+00 ) then call evnnot ( break, coef, l, k, brknew, lnew, coefg ) return end if ! ! For I = 2,..., LNEW, construct BRKNEW(I) = A + G^(-1)(STEPI), ! with STEPI = ( I - 1 ) * STEP. ! ! This requires inversion of the piecewise linear function G. ! ! For this, J is found so that ! ! G(BREAK(J)) <= STEPI <= G(BREAK(J+1)) ! ! and then ! ! BRKNEW(I) = BREAK(J) + ( STEPI - G(BREAK(J)) ) / DG(BREAK(J) ). ! ! The midpoint is chosen if DG(BREAK(J)) = 0. ! j = 1 do i = 2, lnew stepi = real ( i - 1, kind = rk ) * step do if ( j == l ) then exit end if if ( stepi <= coefg(1,j+1) ) then exit end if j = j + 1 end do if ( coefg(2,j) /= 0.0D+00 ) then brknew(i) = break(j) + ( stepi - coefg(1,j) ) / coefg(2,j) else brknew(i) = ( break(j) + break(j+1) ) / 2.0D+00 end if end do return end function ppvalu ( break, coef, l, k, x, jderiv ) !*****************************************************************************80 ! !! PPVALU evaluates a piecewise polynomial function or its derivative. ! ! Discussion: ! ! PPVALU calculates the value at X of the JDERIV-th derivative of ! the piecewise polynomial function F from its piecewise ! polynomial representation. ! ! The interval index I, appropriate for X, is found through a ! call to INTERV. The formula for the JDERIV-th derivative ! of F is then evaluated by nested multiplication. ! ! The J-th derivative of F is given by: ! (d^J) F(X) = ! COEF(J+1,I) + H * ( ! COEF(J+2,I) + H * ( ! ... ! COEF(K-1,I) + H * ( ! COEF(K, I) / (K-J-1) ) / (K-J-2) ... ) / 2 ) / 1 ! with ! H = X - BREAK(I) ! and ! I = max ( 1, max ( J, BREAK(J) <= X, 1 <= J <= L ) ). ! ! Modified: ! ! 16 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BREAK(L+1), real COEF(*), integer L, the ! piecewise polynomial representation of the function F to be evaluated. ! ! Input, integer K, the order of the polynomial pieces that ! make up the function F. The usual value for K is 4, signifying a ! piecewise cubic polynomial. ! ! Input, real ( kind = rk ) X, the point at which to evaluate F or ! of its derivatives. ! ! Input, integer JDERIV, the order of the derivative to be ! evaluated. If JDERIV is 0, then F itself is evaluated, ! which is actually the most common case. It is assumed ! that JDERIV is zero or positive. ! ! Output, real ( kind = rk ) PPVALU, the value of the JDERIV-th ! derivative of F at X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer l real ( kind = rk ) break(l+1) real ( kind = rk ) coef(k,l) real ( kind = rk ) fmmjdr real ( kind = rk ) h integer i integer jderiv integer m integer ndummy real ( kind = rk ) ppvalu real ( kind = rk ) value real ( kind = rk ) x value = 0.0D+00 fmmjdr = k - jderiv ! ! Derivatives of order K or higher are identically zero. ! if ( k <= jderiv ) then ppvalu = value return end if ! ! Find the index I of the largest breakpoint to the left of X. ! call interv ( break, l+1, x, i, ndummy ) ! ! Evaluate the JDERIV-th derivative of the I-th polynomial piece at X. ! h = x - break(i) m = k do value = ( value / fmmjdr ) * h + coef(m,i) m = m - 1 fmmjdr = fmmjdr - 1.0D+00 if ( fmmjdr <= 0.0D+00 ) then exit end if end do ppvalu = value return end subroutine putit ( t, kpm, left, scrtch, dbiatx, q, nrow, b ) !*****************************************************************************80 ! !! PUTIT puts together one block of the collocation equation system. ! ! Method: ! ! The K collocation equations for the interval ( T(LEFT), T(LEFT+1) ) ! are constructed with the aid of the subroutine DIFEQU( 2, ., . ) ! and interspersed (in order) with the side conditions, if any, in ! this interval, using DIFEQU ( 3, ., . ) for the information. ! ! The block Q has KPM columns, corresponding to the KPM B-splines of order ! KPM which have the interval ( T(LEFT), T(LEFT+1) ) in their support. ! ! The block's diagonal is part of the diagonal of the total system. ! ! The first equation in this block not overlapped by the preceding block ! is therefore equation LOWROW, with LOWROW = number of side conditions ! in preceding intervals (or blocks). ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) T(LEFT+KPM), the knot sequence. ! ! Input, integer KPM, the order of the spline. ! ! Input, integer LEFT, indicates the interval of interest, ! that is, the interval ( T(LEFT), T(LEFT+1) ). ! ! Workspace, real ( kind = rk ) SCRTCH(KPM,KPM). ! ! Workspace, real ( kind = rk ) DBIATX(KPM,M+1), derivatives of B-splines, ! with DBIATX(J,I+1) containing the I-th derivative of the J-th B-spline ! of interest. ! ! Output, real ( kind = rk ) Q(NROW,KPM), the block. ! ! Input, integer NROW, number of rows in block to be !! put together. ! ! Output, real ( kind = rk ) B(NROW), the corresponding piece of ! the right hand side. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer kpm integer left integer nrow real ( kind = rk ) b(nrow) real ( kind = rk ) dbiatx(kpm,*) real ( kind = rk ) dx integer irow integer iside integer itermx integer j integer k integer ll integer lowrow integer m integer mode integer mp1 real ( kind = rk ) q(nrow,kpm) real ( kind = rk ) rho real ( kind = rk ) scrtch(kpm,kpm) real ( kind = rk ) t(left+kpm) real ( kind = rk ) v(20) real ( kind = rk ) xm real ( kind = rk ) xside real ( kind = rk ) xx save / other / save / side / common / other / itermx, k, rho(19) common / side / m, iside, xside(10) mp1 = m + 1 q(1:nrow,1:kpm) = 0.0D+00 xm = ( t(left+1) + t(left) ) / 2.0D+00 dx = ( t(left+1) - t(left) ) / 2.0D+00 ll = 1 lowrow = iside do irow = lowrow, nrow if ( k < ll ) then go to 20 end if mode = 2 ! ! Next collocation point: ! xx = xm + dx * rho(ll) ll = ll + 1 ! ! The corresponding collocation equation is next unless the next side ! condition occurs at a point at, or to the left of, the next ! collocation point. ! if ( m < iside ) then go to 30 end if if ( xx < xside(iside) ) then go to 30 end if ll = ll - 1 20 continue mode = 3 xx = xside(iside) 30 continue call difequ ( mode, xx, v ) ! ! The next equation, a collocation equation (MODE=2) or a side ! condition (MODE=3), reads ! ! (*) (V(M+1)*D^M+V(M)*D^(M-1) +...+ V(1)*D^0)F(XX) = V(M+2) ! ! in terms of the information supplied by DIFEQU. ! ! The corresponding equation for the B-spline coefficients of F therefore ! has the left side of (*), evaluated at each of the KPM B-splines having ! XX in their support, as its KPM possibly nonzero coefficients. ! call bsplvd ( t, kpm, xx, left, scrtch, dbiatx, mp1 ) do j = 1, kpm q(irow,j) = dot_product ( dbiatx(j,1:mp1), v(1:mp1) ) end do b(irow) = v(m+2) end do return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is an array of double precision real values. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,g16.8)' ) i, a(i) end do return end function round ( x, size ) !*****************************************************************************80 ! !! ROUND is called to add some noise to data. ! ! Discussion: ! ! This function simply adds plus or minus a perturbation value ! to the input data. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input. real ( kind = rk ) X, the value to be perturbed. ! ! Input, real ( kind = rk ) SIZE, the size of the perturbation. ! ! Output, real ( kind = rk ) ROUND, the perturbed value. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), save :: flip = -1.0D+00 real ( kind = rk ) round real ( kind = rk ) size real ( kind = rk ) x flip = -flip round = x + flip * size return end subroutine sbblok ( bloks, integs, nbloks, ipivot, b, x ) !*****************************************************************************80 ! !! SBBLOK solves a linear system that was factored by FCBLOK. ! ! Discussion: ! ! The routine supervises the solution, by forward and backward ! substitution, of the linear system ! A * x = b ! for X, with the PLU factorization of A already generated in FCBLOK. ! Individual blocks of equations are solved via SUBFOR and SUBBAK. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) BLOKS(*), integer INTEGS(3,NBLOKS), integer ! NBLOKS, integer IPIVOT(*), are as on return from FCBLOK. ! ! Input, real ( kind = rk ) B(*), the right hand side, stored corresponding ! to the storage of the equations. See comments in SLVBLK for details. ! ! Output, real ( kind = rk ) X(*), the solution vector. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbloks real ( kind = rk ) b(*) real ( kind = rk ) bloks(*) integer i integer index integer indexb integer indexx integer integs(3,nbloks) integer ipivot(*) integer j integer last integer nbp1 integer ncol integer nrow real ( kind = rk ) x(*) ! ! Forward substitution: ! index = 1 indexb = 1 indexx = 1 do i = 1, nbloks nrow = integs(1,i) last = integs(3,i) call subfor ( bloks(index), ipivot(indexb), nrow, last, b(indexb), & x(indexx) ) index = nrow * integs(2,i) + index indexb = indexb + nrow indexx = indexx + last end do ! ! Back substitution. ! nbp1 = nbloks + 1 do j = 1, nbloks i = nbp1 - j nrow = integs(1,i) ncol = integs(2,i) last = integs(3,i) index = index - nrow * ncol indexb = indexb - nrow indexx = indexx - last call subbak ( bloks(index), ipivot(indexb), nrow, ncol, last, x(indexx) ) end do return end subroutine setupq ( x, dx, y, npoint, v, qty ) !*****************************************************************************80 ! !! SETUPQ is to be called in SMOOTH. ! ! Discussion: ! ! Put DELX = X(*+1) - X(*) into V(*,4). ! ! Put the three bands of Q' * D into V(*,1:3). ! ! Put the three bands of ( D * Q )' * ( D * Q ) at and above the diagonal ! into V(*,5:7). ! ! Here, Q is the tridiagonal matrix of order ( NPOINT-2, NPOINT ) ! with general row ! 1/DELX(I), -1/DELX(I)-1/DELX(I+1), 1/DELX(I+1) ! and D is the diagonal matrix with general row DX(I). ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) X(NPOINT), the abscissas, assumed to be ! strictly increasing. ! ! Input, real ( kind = rk ) DX(NPOINT), the data uncertainty estimates, ! which are assumed to be positive. ! ! Input, real ( kind = rk ) Y(NPOINT), the corresponding ordinates. ! ! Input, integer NPOINT, the number of data points. ! ! Output, real ( kind = rk ) V(NPOINT,7), contains data needed for ! the smoothing computation. ! ! Output, real ( kind = rk ) QTY(NPOINT), the value of Q' * Y. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer npoint real ( kind = rk ) diff real ( kind = rk ) dx(npoint) integer i real ( kind = rk ) prev real ( kind = rk ) qty(npoint) real ( kind = rk ) v(npoint,7) real ( kind = rk ) x(npoint) real ( kind = rk ) y(npoint) v(1:npoint-1,4) = x(2:npoint) - x(1:npoint-1) v(2:npoint-1,1) = dx(1:npoint-2) / v(1:npoint-2,4) v(npoint,1) = 0.0D+00 v(2:npoint-1,2) = - dx(2:npoint-1) / v(2:npoint-1,4) & - dx(2:npoint-1) / v(1:npoint-2,4) v(2:npoint-1,3) = dx(3:npoint) / v(2:npoint-1,4) v(2:npoint-1,5) = v(2:npoint-1,1)**2 & + v(2:npoint-1,2)**2 & + v(2:npoint-1,3)**2 v(2:npoint-2,6) = v(2:npoint-2,2) * v(3:npoint-1,1) & + v(2:npoint-2,3) * v(3:npoint-1,2) v(npoint-1,6) = 0.0D+00 v(2:npoint-3,7) = v(2:npoint-3,3) * v(4:npoint-1,1) v(npoint-2,7) = 0.0D+00 v(npoint-1,7) = 0.0D+00 ! ! Construct Q' * Y in QTY. ! prev = ( y(2) - y(1) ) / v(1,4) do i = 2, npoint - 1 diff = ( y(i+1) - y(i) ) / v(i,4) qty(i) = diff - prev prev = diff end do return end subroutine shiftb ( ai, ipivot, nrowi, ncoli, last, ai1, nrowi1, ncoli1 ) !*****************************************************************************80 ! !! SHIFTB shifts the rows in the current block. ! ! Discussion: ! ! This routine shifts rows in the current block, AI, which are not used ! as pivot rows, if any, that is, rows IPIVOT(LAST+1) through IPIVOT(NROWI), ! onto the first MMAX = NROW - LAST rows of the next block, AI1, ! with column LAST + J of AI going to column J, ! for J = 1,..., JMAX = NCOLI - LAST. ! ! The remaining columns of these rows of AI1 are zeroed out. ! ! Diagram: ! ! Original situation after Results in a new block I+1 ! LAST = 2 columns have been created and ready to be ! done in FACTRB, assuming no factored by next FACTRB call. ! interchanges of rows. ! ! 1 ! X X 1X X X X X X X X ! 1 ! 0 X 1X X X 0 X X X X ! BLOCK I 1 --------------- ! NROWI=4 0 0 1X X X 0 0 1X X X 0 01 ! NCOLI=5 1 1 1 ! LAST=2 0 0 1X X X 0 0 1X X X 0 01 ! ------------------- 1 1 NEW ! 1X X X X X 1X X X X X1 BLOCK ! 1 1 1 I+1 ! BLOCK I+1 1X X X X X 1X X X X X1 ! NROWI1= 5 1 1 1 ! NCOLI1= 5 1X X X X X 1X X X X X1 ! ------------------- 1-------------1 ! 1 ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) AI(NROWI,NCOLI), the current block. ! ! Input, integer IPIVOT(NROWI), the pivot vector. ! ! Input, integer NROWI, NCOLI, the number of rows and columns ! in block AI. ! ! Input, integer LAST, indicates the last row on which pivoting ! has been carried out. ! ! Input/output, real ( kind = rk ) AI1(NROWI1,NCOLI1), the next block. ! ! Input, integer NROWI1, NCOLI1, the number of rows and columns ! in block AI1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ncoli integer ncoli1 integer nrowi1 integer nrowi real ( kind = rk ) ai(nrowi,ncoli) real ( kind = rk ) ai1(nrowi1,ncoli1) integer ip integer ipivot(nrowi) integer j integer last integer m if ( nrowi - last < 1 ) then return end if if ( ncoli - last < 1 ) then return end if ! ! Put the remainder of block I into AI1. ! do m = 1, nrowi - last ip = ipivot(last+m) do j = 1, ncoli - last ai1(m,j) = ai(ip,last+j) end do end do ! ! Zero out the upper right corner of AI1. ! do j = ncoli + 1 - last, ncoli1 do m = 1, nrowi - last ai1(m,j) = 0.0D+00 end do end do return end subroutine slvblk ( bloks, integs, nbloks, b, ipivot, x, iflag ) !*****************************************************************************80 ! !! SLVBLK solves the almost block diagonal linear system A * x = b. ! ! Discussion: ! ! Such almost block diagonal matrices arise naturally in piecewise ! polynomial interpolation or approximation and in finite element ! methods for two-point boundary value problems. The PLU factorization ! method is implemented here to take advantage of the special structure ! of such systems for savings in computing time and storage requirements. ! ! SLVBLK relies on several auxiliary programs: ! ! FCBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,SCRTCH,IFLAG) ! factors the matrix A. ! ! SBBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,B,X) ! solves the system A*X=B once A is factored. ! ! DTBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,IFLAG,DETSGN,DETLOG) ! computes the determinant of A once it has been factored. ! ! Block structure of A: ! ! The NBLOKS blocks are stored consecutively in the array BLOKS. ! ! The first block has its (1,1)-entry at BLOKS(1), and, if the I-th ! block has its (1,1)-entry at BLOKS(INDEX(I)), then ! ! INDEX(I+1) = INDEX(I) + NROW(I) * NCOL(I). ! ! The blocks are pieced together to give the interesting part of A ! as follows. For I=1,2,..., NBLOKS-1, the (1,1)-entry of the next ! block (the (I+1)st block) corresponds to the (LAST+1,LAST+1)-entry ! of the current I-th block. Recall LAST = INTEGS(3,I) and note that ! this means that ! ! A: every block starts on the diagonal of A. ! ! B: the blocks overlap (usually). the rows of the (I+1)st block ! which are overlapped by the I-th block may be arbitrarily ! defined initially. They are overwritten during elimination. ! ! The right hand side for the equations in the I-th block are stored ! correspondingly as the last entries of a piece of B of length NROW ! (= INTEGS(1,I)) and following immediately in B the corresponding ! piece for the right hand side of the preceding block, with the right ! hand side for the first block starting at B(1). In this, the right ! hand side for an equation need only be specified once on input, ! in the first block in which the equation appears. ! ! Example: ! ! The test driver for this package contains an example, a linear ! system of order 11, whose nonzero entries are indicated in the ! following diagram by their row and column index modulo 10. Next to it ! are the contents of the INTEGS arrray when the matrix is taken to ! be almost block diagonal with NBLOKS = 5, and below it are the five ! blocks. ! ! NROW1 = 3, NCOL1 = 4 ! 11 12 13 14 ! 21 22 23 24 NROW2 = 3, NCOL2 = 3 ! 31 32 33 34 ! LAST1 = 2 43 44 45 ! 53 54 55 NROW3 = 3, NCOL3 = 4 ! LAST2 = 3 66 67 68 69 NROW4 = 3, NCOL4 = 4 ! 76 77 78 79 NROW5 = 4, NCOL5 = 4 ! 86 87 88 89 ! LAST3 = 1 97 98 99 90 ! LAST4 = 1 08 09 00 01 ! 18 19 10 11 ! LAST5 = 4 ! ! Actual input to BLOKS shown by rows of blocks of A. ! The ** items are arbitrary. ! ! 11 12 13 14 / ** ** ** / 66 67 68 69 / ** ** ** ** / ** ** ** ** ! 21 22 23 24 / 43 44 45 / 76 77 78 79 / ** ** ** ** / ** ** ** ** ! 31 32 33 34/ 53 54 55/ 86 87 88 89/ 97 98 99 90/ 08 09 00 01 ! 18 19 10 11 ! ! INDEX = 1 INDEX = 13 INDEX = 22 INDEX = 34 INDEX = 46 ! ! Actual right hand side values with ** for arbitrary values: ! ! B1 B2 B3 ** B4 B5 B6 B7 B8 ** ** B9 ** ** B10 B11 ! ! It would have been more efficient to combine block 3 with block 4. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input/output, real ( kind = rk ) BLOKS(*), a one-dimenional array, ! of length sum ( INTEGS(1,1:NBLOKS) * INTEGS(2,1:NBLOKS) ). ! On input, contains the blocks of the almost block diagonal matrix A. ! The array INTEGS describes the block structure. ! On output, contains correspondingly the PLU factorization ! of A, if IFLAG /= 0. Certain entries in BLOKS are arbitrary, ! where the blocks overlap. ! ! Input, integer INTEGS(3,NBLOKS), description of the block ! structure of A. ! integs(1,I) = number of rows of block I = nrow; ! integs(2,I) = number of colums of block I = ncol; ! integs(3,I) = number of elimination steps in block I = last. ! The linear system is of order n = sum ( integs(3,i), i=1,...,nbloks ), ! but the total number of rows in the blocks is ! nbrows=sum( integs(1,i) ; i = 1,...,nbloks) ! ! Input, integer NBLOKS, the number of blocks. ! ! Input, real ( kind = rk ) B(NBROWS), the right hand side. Certain entries ! are arbitrary, corresponding to rows of the blocks which overlap. See ! the block structure in the example. ! ! Output, integer IPIVOT(NBROWS), the pivoting sequence used. ! ! Output, real ( kind = rk ) X(N), the computed solution, if iflag /= 0. ! ! Output, integer IFLAG. ! = (-1)^(number of interchanges during factorization) if A is invertible; ! = 0 if A is singular. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer nbloks real ( kind = rk ) b(*) real ( kind = rk ) bloks(*) integer iflag integer integs(3,nbloks) integer ipivot(*) real ( kind = rk ) x(*) ! ! In the call to FCBLOK, X is used for temporary storage. ! call fcblok ( bloks, integs, nbloks, ipivot, x, iflag ) if ( iflag == 0 ) then return end if call sbblok ( bloks, integs, nbloks, ipivot, b, x ) return end function smooth ( x, y, dy, npoint, s, v, a ) !*****************************************************************************80 ! !! SMOOTH constructs the cubic smoothing spline to given data. ! ! Discussion: ! ! The data is of the form ! ( X(1:NPOINT), Y(1:NPOINT) ) ! ! The cubic smoothing spline has as small a second derivative as ! possible, while ! S(F) <= S, ! where ! S(F) = sum ( 1 <= I <= NPOINT ) ( ( ( Y(I) - F(X(I)) ) / DY(I) )^2. ! ! Method: ! ! The matrices Q' * D and Q' * D^2 * Q are constructed in SETUPQ from ! X and DY, as is the vector QTY = Q' * Y. ! ! Then, for given P, the vector U is determined in CHOL1D as ! the solution of the linear system ! ( 6 * (1-P) * Q' * D^2 * Q + P * R ) * U = QTY. ! ! From U and this choice of smoothing parameter P, the smoothing spline F ! is obtained in the sense that: ! F(X(.)) = Y - 6 (1-P) D^2 * Q * U, ! (d^2) F(X(.)) = 6 * P * U. ! ! The smoothing parameter P is found, if possible, so that ! SF(P) = S, ! with SF(P) = S(F), where F is the smoothing spline as it depends ! on P. If S = 0, then P = 1. If SF(0) <= S, then P = 0. ! Otherwise, the secant method is used to locate an appropriate P in ! the open interval (0,1). ! ! Specifically, ! P(0) = 0, P(1) = ( S - SF(0) ) / DSF ! with ! DSF = -24 * U' * R * U ! a good approximation to ! D(SF(0)) = DSF + 60 * (D*Q*U)' * (D*Q*U), ! and U as obtained for P = 0. ! ! After that, for N = 1, 2,... until SF(P(N)) <= 1.01 * S, do: ! determine P(N+1) as the point at which the secant to SF at the ! points P(N) and P(N-1) takes on the value S. ! ! If 1 <= P(N+1), choose instead P(N+1) as the point at which ! the parabola SF(P(N))*((1-.)/(1-P(N)))^2 takes on the value S. ! ! Note that, in exact arithmetic, it is always the case that ! P(N+1) < P(N), ! hence ! SF(P(N+1)) < SF(P(N)). ! ! Therefore, also stop the iteration, with final P = 1, in case ! SF(P(N)) <= SF(P(N+1)). ! ! Modified: ! ! 16 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) X(NPOINT), the abscissas, assumed to be ! strictly increasing. ! ! Input, real ( kind = rk ) Y(NPOINT), the corresponding ordinates. ! ! Input, real ( kind = rk ) DY(NPOINT), the data uncertainty estimates, ! which are assumed to be positive. ! ! Input, integer NPOINT, the number of data points. ! ! Input, real ( kind = rk ) S, an upper bound on the discrete weighted mean ! square distance of the approximation F from the data. ! ! Workspace, real ( kind = rk ) V(NPOINT,7). ! ! Workspace, real ( kind = rk ) A(NPOINT,4). ! ! Output, real ( kind = rk ) A(NPOINT,4). ! A(*,1).....contains the sequence of smoothed ordinates. ! A(I,J) = d^(J-1) F(X(I)), for J = 2:4, I = 1:NPOINT-1. ! That is, the first three derivatives of the smoothing spline F at the ! left end of each of the data intervals. Note that A would have to ! be transposed before it could be used in PPVALU. ! ! Output, real ( kind = rk ) SMOOTH, the value of the smoothing parameter. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer npoint real ( kind = rk ) a(npoint,4) real ( kind = rk ) change real ( kind = rk ) dy(npoint) integer i real ( kind = rk ) oosf real ( kind = rk ) ooss real ( kind = rk ) p real ( kind = rk ) prevq real ( kind = rk ) prevsf real ( kind = rk ) q real ( kind = rk ) s real ( kind = rk ) sfq real ( kind = rk ) smooth real ( kind = rk ) utru real ( kind = rk ) v(npoint,7) real ( kind = rk ) x(npoint) real ( kind = rk ) y(npoint) call setupq ( x, dy, y, npoint, v, a(1,4) ) if ( s <= 0.0D+00 ) then p = 1.0D+00 call chol1d ( p, v, a(1,4), npoint, 1, a(1,3), a(1,1) ) sfq = 0.0D+00 else p = 0.0D+00 call chol1d ( p, v, a(1,4), npoint, 1, a(1,3), a(1,1) ) sfq = 36.0D+00 * dot_product ( a(1:npoint,1)**2, dy(1:npoint)**2 ) if ( s < sfq ) then utru = 0.0D+00 do i = 2, npoint utru = utru + v(i-1,4) * ( a(i-1,3) * ( a(i-1,3) + a(i,3) ) & + a(i,3)**2 ) end do ooss = 1.0D+00 / sqrt ( s ) oosf = 1.0D+00 / sqrt ( sfq ) q = - ( oosf - ooss ) * sfq / ( 6.0D+00 * utru * oosf ) ! ! Secant iteration for the determination of P starts here. ! prevq = 0.0D+00 prevsf = oosf do call chol1d ( q / ( 1.0D+00 + q ), v, a(1,4), npoint, 1, & a(1,3), a(1,1) ) sfq = 36.0D+00 * dot_product ( a(1:npoint,1)**2, dy(1:npoint)**2 ) & / ( 1.0D+00 + q )**2 if ( abs ( sfq - s ) <= 0.01D+00 * s ) then exit end if oosf = 1.0D+00 / sqrt ( sfq ) change = ( q - prevq ) / ( oosf - prevsf ) * ( oosf - ooss ) prevq = q q = q - change prevsf = oosf end do p = q / ( 1.0D+00 + q ) end if end if ! ! Correct value of P has been found. ! Compute polynomial coefficients from Q * U in A(.,1). ! smooth = sfq a(1:npoint,1) = y(1:npoint) - 6.0D+00 * ( 1.0D+00 - p ) & * dy(1:npoint)**2 * a(1:npoint,1) a(1:npoint,3) = a(1:npoint,3) * 6.0D+00 * p do i = 1, npoint - 1 a(i,4) = ( a(i+1,3) - a(i,3) ) / v(i,4) a(i,2) = ( a(i+1,1) - a(i,1) ) / v(i,4) & - ( a(i,3) + a(i,4) / 3.0D+00 * v(i,4) ) / 2.0D+00 * v(i,4) end do return end subroutine spli2d ( tau, gtau, t, n, k, m, work, q, bcoef, iflag ) !*****************************************************************************80 ! !! SPLI2D produces a interpolatory tensor product spline. ! ! Discussion: ! ! SPLI2D is an extended version of SPLINT. ! ! SPLI2D produces the B-spline coefficients BCOEF(J,.) of the ! spline of order K with knots T(1:N+K), which takes on ! the value GTAU(I,J) at TAU(I), I=1,..., N, J=1,...,M. ! ! The I-th equation of the linear system ! A * BCOEF = B ! for the B-spline coefficients of the interpolant enforces ! interpolation at TAU(I), I=1,...,N. Hence, B(I) = GTAU(I), ! for all I, and A is a band matrix with 2*K-1 bands, if it is ! invertible. ! ! The matrix A is generated row by row and stored, diagonal by ! diagonal, in the rows of the array Q, with the main diagonal ! going into row K. ! ! The banded system is then solved by a call to BANFAC, which ! constructs the triangular factorization for A and stores it ! again in Q, followed by a call to BANSLV, which then obtains ! the solution BCOEF by substitution. ! ! The linear system to be solved is theoretically invertible if ! and only if ! T(I) < TAU(I) < TAU(I+K), for all I. ! Violation of this condition is certain to lead to IFLAG = 2. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) TAU(N), contains the data point abscissas. ! TAU must be strictly increasing ! ! Input, real ( kind = rk ) GTAU(N,M), contains the data point ordinates. ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. ! ! Input, integer N, the number of data points and the ! dimension of the spline space SPLINE(K,T) ! ! Input, integer K, the order of the spline. ! ! Input, integer M, the number of data sets. ! ! Work space, real ( kind = rk ) WORK(N). ! ! Output, real ( kind = rk ) Q(2*K-1)*N, the triangular ! factorization of the coefficient matrix of the linear ! system for the B-spline coefficients of the spline interpolant. ! The B-spline coefficients for the interpolant of an additional ! data set ( TAU(I), HTAU(I) ), I=1,...,N with the same data ! abscissae can be obtained without going through all the ! calculations in this routine, simply by loading HTAU into ! BCOEF and then using the statement ! CALL BANSLV ( Q, 2*K-1, N, K-1, K-1, BCOEF ) ! ! Output, real ( kind = rk ) BCOEF(N), the B-spline coefficients of ! the interpolant. ! ! Output, integer IFLAG, error indicator. ! 1, no error. ! 2, an error occurred, which may have been caused by ! singularity of the linear system. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) bcoef(m,n) real ( kind = rk ) gtau(n,m) integer i integer iflag integer ilp1mx integer j integer jj integer k integer left real ( kind = rk ) q((2*k-1)*n) real ( kind = rk ) t(n+k) real ( kind = rk ) tau(n) real ( kind = rk ) taui real ( kind = rk ) work(n) left = k q(1:(2*k-1)*n) = 0.0D+00 ! ! Construct the N interpolation equations. ! do i = 1, n taui = tau(i) ilp1mx = min ( i + k, n + 1 ) ! ! Find the index LEFT in the closed interval (I,I+K-1) such that: ! ! T(LEFT) < = TAU(I) < T(LEFT+1) ! ! The matrix will be singular if this is not possible. ! left = max ( left, i ) if ( taui < t(left) ) then iflag = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLI2D - Fatal error!' write ( *, '(a)' ) ' The TAU array is not strictly increasing.' stop 1 end if do while ( t(left+1) <= taui ) left = left + 1 if ( left < ilp1mx ) then cycle end if left = left - 1 if ( t(left+1) < taui ) then iflag = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLI2D - Fatal error!' write ( *, '(a)' ) ' The TAU array is not strictly increasing.' stop 1 end if exit end do ! ! The I-th equation enforces interpolation at TAUI, hence ! ! A(I,J) = B(J,K,T)(TAUI), for all J. ! ! Only the K entries with J = LEFT-K+1, ..., LEFT actually might be ! nonzero. These K numbers are returned, in WORK (used for ! temporary storage here), by the following call: ! call bsplvb ( t, k, 1, taui, left, work ) ! ! We therefore want ! WORK(J) = B(LEFT-K+J)(TAUI) ! to go into ! A(I,LEFT-K+J), ! ! that is, into Q(I-(LEFT+J)+2*K,(LEFT+J)-K) since ! A(I+J,J) is to go into Q(I+K,J), for all I, J, if we consider Q ! as a two-dimensional array, with 2*K-1 rows. See comments in ! BANFAC. ! ! In the present program, we treat Q as an equivalent one-dimensional ! array, because of fortran restrictions on dimension statements. ! ! We therefore want WORK(J) to go into the entry of Q with index: ! I -(LEFT+J)+2*K + ((LEFT+J)-K-1)*(2*K-1) ! = I-LEFT+1+(LEFT -K)*(2*K-1) + (2*K-2)*J ! jj = i - left + 1 + ( left - k ) * ( k + k - 1 ) do j = 1, k jj = jj + k + k - 2 q(jj) = work(j) end do end do ! ! Factor A, stored again in Q. ! call banfac ( q, k+k-1, n, k-1, k-1, iflag ) if ( iflag == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLI2D - Fatal error!' write ( *, '(a)' ) ' BANFAC reports that the matrix is singular.' stop 1 end if ! ! Solve ! A * BCOEF = GTAU ! by back substitution. ! do j = 1, m work(1:n) = gtau(1:n,j) call banslv ( q, k+k-1, n, k-1, k-1, work ) bcoef(j,1:n) = work(1:n) end do return end subroutine splint ( tau, gtau, t, n, k, q, bcoef, iflag ) !*****************************************************************************80 ! !! SPLINT produces the B-spline coefficients BCOEF of an interpolating spline. ! ! Discussion: ! ! The spline is of order K with knots T(1:N+K), and takes on the ! value GTAU(I) at TAU(I), for I = 1 to N. ! ! The I-th equation of the linear system ! A * BCOEF = B ! for the B-spline coefficients of the interpolant enforces interpolation ! at TAU(1:N). ! ! Hence, B(I) = GTAU(I), for all I, and A is a band matrix with 2*K-1 ! bands, if it is invertible. ! ! The matrix A is generated row by row and stored, diagonal by diagonal, ! in the rows of the array Q, with the main diagonal going ! into row K. See comments in the program. ! ! The banded system is then solved by a call to BANFAC, which ! constructs the triangular factorization for A and stores it again in ! Q, followed by a call to BANSLV, which then obtains the solution ! BCOEF by substitution. ! ! BANFAC does no pivoting, since the total positivity of the matrix ! A makes this unnecessary. ! ! The linear system to be solved is (theoretically) invertible if ! and only if ! T(I) < TAU(I) < TAU(I+K), for all I. ! Violation of this condition is certain to lead to IFLAG = 2. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) TAU(N), the data point abscissas. The entries in ! TAU should be strictly increasing. ! ! Input, real ( kind = rk ) GTAU(N), the data ordinates. ! ! Input, real ( kind = rk ) T(N+K), the knot sequence. ! ! Input, integer N, the number of data points. ! ! Input, integer K, the order of the spline. ! ! Output, real ( kind = rk ) Q((2*K-1)*N), the triangular factorization ! of the coefficient matrix of the linear system for the B-coefficients ! of the spline interpolant. The B-coefficients for the interpolant ! of an additional data set can be obtained without going through all ! the calculations in this routine, simply by loading HTAU into BCOEF ! and then executing the call: ! call banslv ( q, 2*k-1, n, k-1, k-1, bcoef ) ! ! Output, real ( kind = rk ) BCOEF(N), the B-spline coefficients of ! the interpolant. ! ! Output, integer IFLAG, error flag. ! 1, = success. ! 2, = failure. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) bcoef(n) real ( kind = rk ) gtau(n) integer i integer iflag integer ilp1mx integer j integer jj integer k integer kpkm2 integer left real ( kind = rk ) q((2*k-1)*n) real ( kind = rk ) t(n+k) real ( kind = rk ) tau(n) real ( kind = rk ) taui kpkm2 = 2 * ( k - 1 ) left = k q(1:(2*k-1)*n) = 0.0D+00 ! ! Loop over I to construct the N interpolation equations. ! do i = 1, n taui = tau(i) ilp1mx = min ( i + k, n + 1 ) ! ! Find LEFT in the closed interval (I,I+K-1) such that ! ! T(LEFT) <= TAU(I) < T(LEFT+1) ! ! The matrix is singular if this is not possible. ! left = max ( left, i ) if ( taui < t(left) ) then iflag = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINT - Fatal Error!' write ( *, '(a)' ) ' The linear system is not invertible!' return end if do while ( t(left+1) <= taui ) left = left + 1 if ( left < ilp1mx ) then cycle end if left = left - 1 if ( t(left+1) < taui ) then iflag = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINT - Fatal Error!' write ( *, '(a)' ) ' The linear system is not invertible!' return end if exit end do ! ! The I-th equation enforces interpolation at TAUI, hence for all J, ! A(I,J) = B(J,K,T)(TAUI). ! ! Only the K entries with J = LEFT-K+1,...,LEFT actually might be nonzero. ! ! These K numbers are returned, in BCOEF (used for temporary storage here), ! by the following. ! call bsplvb ( t, k, 1, taui, left, bcoef ) ! ! We therefore want BCOEF(J) = B(LEFT-K+J)(TAUI) to go into ! A(I,LEFT-K+J), that is, into Q(I-(LEFT+J)+2*K,(LEFT+J)-K) since ! A(I+J,J) is to go into Q(I+K,J), for all I, J, if we consider Q ! as a two-dimensional array, with 2*K-1 rows. See comments in ! BANFAC. ! ! In the present program, we treat Q as an equivalent ! one-dimensional array, because of fortran restrictions on ! dimension statements. ! ! We therefore want BCOEF(J) to go into the entry of Q with index: ! ! I -(LEFT+J)+2*K + ((LEFT+J)-K-1)*(2*K-1) ! = I-LEFT+1+(LEFT -K)*(2*K-1) + (2*K-2)*J ! jj = i - left + 1 + ( left - k ) * ( k + k - 1 ) do j = 1, k jj = jj + kpkm2 q(jj) = bcoef(j) end do end do ! ! Obtain factorization of A, stored again in Q. ! call banfac ( q, k+k-1, n, k-1, k-1, iflag ) if ( iflag == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLINT - Fatal Error!' write ( *, '(a)' ) ' The linear system is not invertible!' return end if ! ! Solve ! ! A * BCOEF = GTAU ! ! by back substitution. ! bcoef(1:n) = gtau(1:n) call banslv ( q, k+k-1, n, k-1, k-1, bcoef ) return end subroutine splopt ( tau, n, k, scrtch, t, iflag ) !*****************************************************************************80 ! !! SPLOPT computes the knots for an optimal recovery scheme. ! ! Discussion: ! ! The optimal recovery scheme is of order K for data at TAU(1:N). ! ! The interior knots T(K+1:N) are determined by Newton's method in ! such a way that the signum function which changes sign at ! T(K+1:N) and nowhere else in ( TAU(1), TAU(N) ) is ! orthogonal to the spline space SPLINE ( K, TAU ) on that interval. ! ! Let XI(J) be the current guess for T(K+J), J=1,...,N-K. Then ! the next Newton iterate is of the form ! XI(J) + (-1)^(N-K-J)*X(J), J=1,...,N-K, ! with X the solution of the linear system ! C * X = D. ! ! Here, for all J, ! C(I,J) = B(I)(XI(J)), ! with B(I) the I-th B-spline of order K for the knot sequence TAU, ! for all I, and D is the vector given, for each I, by ! D(I) = sum ( -A(J), J=I,...,N ) * ( TAU(I+K) - TAU(I) ) / K, ! with, for I = 1 to N-1: ! A(I) = sum ( (-1)^(N-K-J)*B(I,K+1,TAU)(XI(J)), J=1,...,N-K ) ! and ! A(N) = -0.5. ! ! See Chapter XIII of text and references there for a derivation. ! ! The first guess for T(K+J) is sum ( TAU(J+1:J+K-1) ) / ( K - 1 ). ! ! The iteration terminates if max ( abs ( X(J) ) ) < TOL, with ! TOL = TOLRTE * ( TAU(N) - TAU(1) ) / ( N - K ), ! or else after NEWTMX iterations. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) TAU(N), the interpolation points. ! assumed to be nondecreasing, with TAU(I) < TAU(I+K), for all I. ! ! Input, integer N, the number of data points. ! ! Input, integer K, the order of the optimal recovery scheme ! to be used. ! ! Workspace, real ( kind = rk ) SCRTCH((N-K)*(2*K+3)+5*K+3). The various ! contents are specified in the text below. ! ! Output, real ( kind = rk ) T(N+K), the optimal knots ready for ! use in optimal recovery. Specifically, T(1:K) = TAU(1), ! T(N+1:N+K) = TAU(N), while the N - K interior knots T(K+1:N) ! are calculated. ! ! Output, integer IFLAG, error indicator. ! = 1, success. T contains the optimal knots. ! = 2, failure. K < 3 or N < K or the linear system was singular. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer k integer n real ( kind = rk ) del real ( kind = rk ) delmax real ( kind = rk ) floatk integer i integer id integer iflag integer index integer j integer kp1 integer kpkm1 integer kpn integer left integer leftmk integer lenw integer ll integer llmax integer llmin integer na integer nb integer nc integer nd integer, parameter :: newtmx = 10 integer newton integer nmk integer nx real ( kind = rk ) scrtch((n-k)*(2*k+3)+5*k+3) real ( kind = rk ) t(n+k) real ( kind = rk ) tau(n) real ( kind = rk ) sign real ( kind = rk ) signst real ( kind = rk ) tol real ( kind = rk ), parameter :: tolrte = 0.000001D+00 real ( kind = rk ) xij nmk = n - k if ( n < k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLOPT - Fatal error!' write ( *, '(a)' ) ' N < K.' iflag = 2 return end if if ( n == k ) then t(1:k) = tau(1) t(n+1:n+k) = tau(n) return end if if ( k <= 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLOPT - Fatal error!' write ( *, '(a)' ) ' K < 2.' iflag = 2 stop 1 end if floatk = k kp1 = k + 1 kpkm1 = k + k - 1 kpn = k + n signst = -1.0D+00 if ( ( nmk / 2 ) * 2 < nmk ) then signst = 1.0D+00 end if ! ! SCRTCH(I) = TAU-EXTENDED(I), I=1,...,N+K+K ! nx = n + k + k + 1 ! ! SCRTCH(I+NX) = XI(I), I=0,...,N-K+1 ! na = nx + nmk + 1 ! ! SCRTCH(I+NA) = - A(I), I=1,...,N ! nd = na + n ! ! SCRTCH(I+ND) = X(I) or D(I), I=1,...,N-K ! nb = nd + nmk ! ! SCRTCH(I+NB) = BIATX(I), I=1,...,K+1 ! nc = nb + kp1 ! ! SCRTCH(I+(J-1)*(2K-1)+NC) = W(I,J) = C(I-K+J,J), I=J-K,...,J+K, ! J=1,...,N-K. ! lenw = kpkm1 * nmk ! ! Extend TAU to a knot sequence and store in SCRTCH. ! scrtch(1:k) = tau(1) scrtch(k+1:k+n) = tau(1:n) scrtch(kpn+1:kpn+k) = tau(n) ! ! First guess for SCRTCH (.+NX) = XI. ! scrtch(nx) = tau(1) scrtch(nmk+1+nx) = tau(n) do j = 1, nmk scrtch(j+nx) = sum ( tau(j+1:j+k-1) ) / real ( k - 1, kind = rk ) end do ! ! Last entry of SCRTCH (.+NA) = -A is always ... ! scrtch(n+na) = 0.5D+00 ! ! Start the Newton iteration. ! newton = 1 tol = tolrte * ( tau(n) - tau(1) ) / real ( nmk, kind = rk ) ! ! Start the Newton step. ! Compute the 2*K-1 bands of the matrix C and store in SCRTCH(.+NC), ! and compute the vector SCRTCH(.+NA) = -A. ! do newton = 1, newtmx scrtch(nc+1:nc+lenw) = 0.0D+00 scrtch(na+1:na+n-1) = 0.0D+00 sign = signst left = kp1 do j = 1, nmk xij = scrtch(j+nx) do if ( xij < scrtch(left+1) ) then exit end if left = left + 1 if ( kpn <= left ) then left = left - 1 exit end if end do call bsplvb ( scrtch, k, 1, xij, left, scrtch(1+nb) ) ! ! The TAU sequence in SCRTCH is preceded by K additional knots. ! ! Therefore, SCRTCH(LL+NB) now contains B(LEFT-2K+LL)(XIJ) ! which is destined for C(LEFT-2K+LL,J), and therefore for ! ! W(LEFT-K-J+LL,J)= SCRTCH(LEFT-K-J+LL+(J-1)*KPKM1 + NC) ! ! since we store the 2*K-1 bands of C in the 2*K-1 rows of ! the work array W, and W in turn is stored in SCRTCH, ! with W(1,1) = SCRTCH(1+NC). ! ! Also, C being of order N - K, we would want ! 1 <= LEFT-2K+LL <= N - K or ! LLMIN=2K-LEFT <= LL <= N-LEFT+K = LLMAX. ! leftmk = left - k index = leftmk - j + ( j - 1 ) * kpkm1 + nc llmin = max ( 1, k - leftmk ) llmax = min ( k, n - leftmk ) do ll = llmin, llmax scrtch(ll+index) = scrtch(ll+nb) end do call bsplvb ( scrtch, kp1, 2, xij, left, scrtch(1+nb) ) id = max ( 0, leftmk - kp1 ) llmin = 1 - min ( 0, leftmk - kp1 ) do ll = llmin, kp1 id = id + 1 scrtch(id+na) = scrtch(id+na) - sign * scrtch(ll+nb) end do sign = - sign end do call banfac ( scrtch(1+nc), kpkm1, nmk, k-1, k-1, iflag ) if ( iflag == 2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLOPT - Fatal error!' write ( *, '(a)' ) ' Matrix C is not invertible.' stop 1 end if ! ! Compute SCRTCH(.+ND) = D from SCRTCH(.+NA) = -A. ! do i = n, 2, -1 scrtch(i-1+na) = scrtch(i-1+na) + scrtch(i+na) end do do i = 1, nmk scrtch(i+nd) = scrtch(i+na) * ( tau(i+k) - tau(i) ) / floatk end do ! ! Compute SCRTCH(.+ND)= X. ! call banslv ( scrtch(1+nc), kpkm1, nmk, k-1, k-1, scrtch(1+nd) ) ! ! Compute SCRTCH(.+ND) = change in XI. Modify, if necessary, to ! prevent new XI from moving more than 1/3 of the way to its ! neighbors. Then add to XI to obtain new XI in SCRTCH(.+NX). ! delmax = 0.0D+00 sign = signst do i = 1, nmk del = sign * scrtch(i+nd) delmax = max ( delmax, abs ( del ) ) if ( 0.0D+00 < del ) then del = min ( del, ( scrtch(i+1+nx) - scrtch(i+nx) ) / 3.0D+00 ) else del = max ( del, ( scrtch(i-1+nx) - scrtch(i+nx) ) / 3.0D+00 ) end if sign = - sign scrtch(i+nx) = scrtch(i+nx) + del end do ! ! Call it a day in case change in XI was small enough or too many ! steps were taken. ! if ( delmax < tol ) then exit end if end do if ( tol <= delmax ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SPLOPT - Warning!' write ( *, '(a)' ) ' The Newton iteration did not converge.' end if t(1:k) = tau(1) t(k+1:n) = scrtch(nx+1:nx+n-k) t(n+1:n+k) = tau(n) return end subroutine subbak ( w, ipivot, nrow, ncol, last, x ) !*****************************************************************************80 ! !! SUBBAK carries out back substitution for the current block. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) W(NROW,NCOL), integer IPIVOT(NROW), integer ! NROW, integer NCOL, integer LAST, are as on return from FACTRB. ! ! Input/output, real ( kind = rk ) X(NCOL). ! On input, the right hand side for the equations in this block after ! back substitution has been carried out up to, but not including, ! equation IPIVOT(LAST). This means that X(1:LAST) contains the right hand ! sides of equation IPIVOT(1:LAST) as modified during elimination, ! while X(LAST+1:NCOL) is already a component of the solution vector. ! On output, the components of the solution corresponding to the present ! block. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ncol integer nrow integer ip integer ipivot(nrow) integer k integer last real ( kind = rk ) w(nrow,ncol) real ( kind = rk ) x(ncol) do k = last, 1, -1 ip = ipivot(k) x(k) = ( x(k) - dot_product ( w(ip,k+1:ncol), x(k+1:ncol) ) ) / w(ip,k) end do return end subroutine subfor ( w, ipivot, nrow, last, b, x ) !*****************************************************************************80 ! !! SUBFOR carries out the forward pass of substitution for the current block. ! ! Discussion: ! ! The forward pass is the action on the right hand side corresponding to the ! elimination carried out in FACTRB for this block. ! ! At the end, X(1:NROW) contains the right hand side of the transformed ! IPIVOT(1:NROW)-th equation in this block. ! ! Then, since for I=1,...,NROW-LAST, B(NROW+I) is going to be used as ! the right hand side of equation I in the next block (shifted over there ! from this block during factorization), it is set equal to X(LAST+I) here. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) W(NROW,LAST), integer IPIVOT(NROW), ! integer NROW, integer LAST, are as on return from FACTRB. ! ! Output, real ( kind = rk ) B(2*NROW-LAST). On input, B(1:NROW) ! contains the right hand sides for this block. On output, ! B(NROW+1:2*NROW-LAST) contains the appropriately modified right ! hand sides for the next block. ! ! Output, real X(NROW), contains, on output, the appropriately modified ! right hand sides of equations IPIVOT(1:NROW). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer last integer nrow real ( kind = rk ) b(nrow+nrow-last) integer ip integer ipivot(nrow) integer jhi integer k real ( kind = rk ) w(nrow,last) real ( kind = rk ) x(nrow) ip = ipivot(1) x(1) = b(ip) do k = 2, nrow ip = ipivot(k) jhi = min ( k - 1, last ) x(k) = b(ip) - dot_product ( w(ip,1:jhi), x(1:jhi) ) end do ! ! Transfer modified right hand sides of equations IPIVOT(LAST+1:NROW) ! to next block. ! b(nrow+1:2*nrow-last) = x(last+1:nrow) return end subroutine tautsp ( tau, gtau, ntau, gamma, s, break, coef, l, k, iflag ) !*****************************************************************************80 ! !! TAUTSP constructs a cubic spline interpolant to given data. ! ! Discussion: ! ! If 0 < GAMMA, additional knots are introduced where needed to ! make the interpolant more flexible locally. This avoids extraneous ! inflection points typical of cubic spline interpolation at knots to ! rapidly changing data. ! ! Method: ! ! On the I-th interval, (TAU(I), TAU(I+1)), the interpolant is of the ! form: ! (*) F(U(X)) = A + B * U + C * H(U,Z) + D * H(1-U,1-Z), ! with ! U = U(X) = ( X - TAU(I) ) / DTAU(I). ! ! Here, ! Z(I) = ADDG(I+1) / ( ADDG(I) + ADDG(I+1) ) ! but if the denominator vanishes, we set Z(I) = 0.5 ! ! Also, we have ! ADDG(J) = abs ( DDG(J) ), ! DDG(J) = DG(J+1) - DG(J), ! DG(J) = DIVDIF(J) = ( GTAU(J+1) - GTAU(J) ) / DTAU(J) ! and ! H(U,Z) = ALPHA * U^3 ! + ( 1 - ALPHA ) * ( max ( ( ( U - ZETA ) / ( 1 - ZETA ) ), 0 )^3 ! with ! ALPHA(Z) = ( 1 - GAMMA / 3 ) / ZETA ! ZETA(Z) = 1 - GAMMA * min ( ( 1 - Z ), 1/3 ) ! ! Thus, for 1/3 <= Z <= 2/3, F is just a cubic polynomial on ! the interval I. Otherwise, it has one additional knot, at ! TAU(I) + ZETA * DTAU(I). ! ! As Z approaches 1, H(.,Z) has an increasingly sharp bend near 1, ! thus allowing F to turn rapidly near the additional knot. ! ! In terms of F(J) = GTAU(J) and FSECND(J) = second derivative of F ! at TAU(J), the coefficients for (*) are given as: ! A = F(I) - D ! B = ( F(I+1) - F(I) ) - ( C - D ) ! C = FSECND(I+1) * DTAU(I)^2 / HSECND(1,Z) ! D = FSECND(I) * DTAU(I)^2 / HSECND(1,1-Z) ! ! Hence these can be computed once FSECND(1:NTAU) is fixed. ! ! F is automatically continuous and has a continuous second derivative ! except when Z=0 or 1 for some I. We determine FSECND from ! the requirement that the first derivative of F be continuous. ! ! In addition, we require that the third derivative be continuous ! across TAU(2) and across TAU(NTAU-1). This leads to a strictly ! diagonally dominant tridiagonal linear system for the FSECND(I) ! which we solve by Gauss elimination without pivoting. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Input, real ( kind = rk ) TAU(NTAU), the sequence of data points. ! TAU must be strictly increasing. ! ! Input, real ( kind = rk ) GTAU(NTAU), the corresponding sequence of ! function values. ! ! Input, integer NTAU, the number of data points. ! NTAU must be at least 4. ! ! Input, real ( kind = rk ) GAMMA, indicates whether additional flexibility ! is desired. ! GAMMA = 0.0, no additional knots; ! GAMMA in (0.0,3.0), under certain conditions on the given data at ! points I-1, I, I+1, and I+2, a knot is added in the I-th interval, ! for I = 2,...,NTAU-2. See description of method. The interpolant ! gets rounded with increasing gamma. A value of 2.5 for GAMMA is typical. ! GAMMA in (3.0,6.0), same, except that knots might also be added in ! intervals in which an inflection point would be permitted. A value ! of 5.5 for GAMMA is typical. ! ! Output, real ( kind = rk ) BREAK(L), real ( kind = rk ) COEF(K,L), ! integer L, integer K, give the piecewise polynomial ! representation of the interpolant. Specifically, ! for BREAK(i) <= X <= BREAK(I+1), the interpolant has the form: ! F(X) = COEF(1,I) + DX * ( ! COEF(2,I) + (DX/2) * ( ! COEF(3,I) + (DX/3) * ! COEF(4,I) ) ) ! with DX = X - BREAK(I) for I = 1,..., L. ! ! Output, integer IFLAG, error flag. ! 1, no error. ! 2, input was incorrect. ! ! Output, real ( kind = rk ) S(NTAU,6). The individual columns of this ! array contain the following quantities mentioned in the write up ! and below. ! S(.,1) = DTAU = TAU(.+1)-TAU; ! S(.,2) = DIAG = diagonal in linear system; ! S(.,3) = U = upper diagonal in linear system; ! S(.,4) = R = right hand side for linear system (initially) ! = FSECND = solution of linear system, namely the second ! derivatives of interpolant at TAU; ! S(.,5) = Z = indicator of additional knots; ! S(.,6) = 1/HSECND(1,X) with X = Z or 1-Z. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ntau real ( kind = rk ) alph real ( kind = rk ) alpha real ( kind = rk ) break(*) real ( kind = rk ) c real ( kind = rk ) coef(4,*) real ( kind = rk ) d real ( kind = rk ) del real ( kind = rk ) denom real ( kind = rk ) divdif real ( kind = rk ) entry real ( kind = rk ) entry3 real ( kind = rk ) factor real ( kind = rk ) factr2 real ( kind = rk ) gam real ( kind = rk ) gamma real ( kind = rk ) gtau(ntau) integer i integer iflag integer k integer l integer method real ( kind = rk ) onemg3 real ( kind = rk ) onemzt real ( kind = rk ) ratio real ( kind = rk ) s(ntau,6) real ( kind = rk ) sixth real ( kind = rk ) tau(ntau) real ( kind = rk ) temp real ( kind = rk ) x real ( kind = rk ) z real ( kind = rk ) zeta real ( kind = rk ) zt2 alph(x) = min ( 1.0D+00, onemg3 / x ) ! ! There must be at least 4 interpolation points. ! if ( ntau < 4 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TAUTSP - Fatal error!' write ( *, '(a)' ) ' Input NTAU must be at least 4.' write ( *, '(a,i8)' ) ' NTAU = ', ntau iflag = 2 stop 1 end if ! ! Construct delta TAU and first and second (divided) differences of data. ! do i = 1, ntau - 1 s(i,1) = tau(i+1) - tau(i) if ( s(i,1) <= 0.0D+00 ) then write ( *, '(a,i3,a,2e15.6,a)' ) & ' Point ', i, ' and the next ', tau(i), tau(i+1), ' are disordered.' iflag = 2 return end if s(i+1,4) = ( gtau(i+1) - gtau(i) ) / s(i,1) end do do i = 2, ntau - 1 s(i,4) = s(i+1,4) - s(i,4) end do ! ! Construct system of equations for second derivatives at TAU. ! ! At each interior data point, there is one continuity equation. ! At the first and last interior data point there is an additional ! equation for a total of NTAU equations in NTAU unknowns. ! i = 2 s(2,2) = s(1,1) / 3.0D+00 sixth = 1.0D+00 / 6.0D+00 method = 2 gam = gamma if ( gam <= 0.0D+00 ) then method = 1 end if if ( 3.0D+00 < gam ) then method = 3 gam = gam - 3.0D+00 end if onemg3 = 1.0D+00 - gam / 3.0D+00 ! ! Loop over I. ! 70 continue ! ! Construct Z(I) and ZETA(I). ! z = 0.5D+00 if ( method == 1 ) then go to 100 end if if ( method == 3 ) then go to 90 end if if ( s(i,4) * s(i+1,4) < 0.0D+00 ) then go to 100 end if 90 continue temp = abs ( s(i+1,4) ) denom = abs ( s(i,4) ) + temp if ( denom /= 0.0D+00 ) then z = temp / denom if ( abs ( z - 0.5D+00 ) <= sixth ) then z = 0.5D+00 end if end if 100 continue s(i,5) = z ! ! Set up part of the I-th equation which depends on the I-th interval. ! if ( z < 0.5D+00 ) then zeta = gam * z onemzt = 1.0D+00 - zeta zt2 = zeta**2 alpha = alph(onemzt) factor = zeta / ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 ) s(i,6) = zeta * factor / 6.0D+00 s(i,2) = s(i,2) + s(i,1) & * ( ( 1.0D+00 - alpha * onemzt ) * factor / 2.0D+00 - s(i,6) ) ! ! If Z = 0 and the previous Z = 1, then D(I) = 0. ! Since then also U(I-1) = L(I+1) = 0, its value does not matter. ! Reset D(I) = 1 to insure nonzero pivot in elimination. ! if ( s(i,2) <= 0.0D+00 ) then s(i,2) = 1.0D+00 end if s(i,3) = s(i,1) / 6.0D+00 else if ( z == 0.5D+00 ) then s(i,2) = s(i,2) + s(i,1) / 3.0D+00 s(i,3) = s(i,1) / 6.0D+00 else if ( 0.5D+00 < z ) then onemzt = gam * ( 1.0D+00 - z ) zeta = 1.0D+00 - onemzt alpha = alph(zeta) factor = onemzt / ( 1.0D+00 - alpha * zeta * ( 1.0D+00 + onemzt ) ) s(i,6) = onemzt * factor / 6.0D+00 s(i,2) = s(i,2) + s(i,1) / 3.0D+00 s(i,3) = s(i,6) * s(i,1) end if if ( 2 < i ) then go to 190 end if s(1,5) = 0.5D+00 ! ! The first two equations enforce continuity of the first and of ! the third derivative across TAU(2). ! s(1,2) = s(1,1) / 6.0D+00 s(1,3) = s(2,2) entry3 = s(2,3) if ( z < 0.5D+00 ) then factr2 = zeta * ( alpha * ( zt2 - 1.0D+00 ) + 1.0D+00 ) & / ( alpha * ( zeta * zt2 - 1.0D+00 ) + 1.0D+00 ) ratio = factr2 * s(2,1) / s(1,2) s(2,2) = factr2 * s(2,1) + s(1,1) s(2,3) = - factr2 * s(1,1) else if ( z == 0.5D+00 ) then ratio = s(2,1) / s(1,2) s(2,2) = s(2,1) + s(1,1) s(2,3) = - s(1,1) else if ( 0.5D+00 < z ) then ratio = s(2,1) / s(1,2) s(2,2) = s(2,1) + s(1,1) s(2,3) = - s(1,1) * 6.0D+00 * alpha * s(2,6) end if ! ! At this point, the first two equations read: ! DIAG(1)*X1+U(1)*X2 + ENTRY3*X3 = R(2) ! -RATIO*DIAG(1)*X1+DIAG(2)*X2 + U(2)*X3 = 0.0 ! Eliminate first unknown from second equation. ! s(2,2) = ratio * s(1,3) + s(2,2) s(2,3) = ratio * entry3 + s(2,3) s(1,4) = s(2,4) s(2,4) = ratio * s(1,4) go to 200 190 continue ! ! The I-th equation enforces continuity of the first derivative ! across TAU(I). It now reads: ! - RATIO * DIAG(I-1) * X(I-1) + DIAG(I) * X(I) + U(I) * X(I+1) = R(I). ! Eliminate (I-1)st unknown from this equation ! s(i,2) = ratio * s(i-1,3) + s(i,2) s(i,4) = ratio * s(i-1,4) + s(i,4) ! ! Set up the part of the next equation which depends on the I-th interval. ! 200 continue if ( z < 0.5D+00 ) then ratio = - s(i,6) * s(i,1) / s(i,2) s(i+1,2) = s(i,1) / 3.0D+00 else if ( z == 0.5D+00 ) then ratio = - ( s(i,1) / 6.0D+00 ) / s(i,2) s(i+1,2) = s(i,1) / 3.0D+00 else if ( 0.5D+00 < z ) then ratio = - ( s(i,1) / 6.0D+00 ) / s(i,2) s(i+1,2) = s(i,1) & * ( ( 1.0D+00 - zeta * alpha ) * factor / 2.0D+00 - s(i,6) ) end if ! ! End of I loop. ! i = i + 1 if ( i < ntau - 1 ) then go to 70 end if s(i,5) = 0.5D+00 ! ! The last two equations enforce continuity of third derivative and ! of first derivative across TAU(NTAU-1). ! entry = ratio * s(i-1,3) + s(i,2) + s(i,1) / 3.0D+00 s(i+1,2) = s(i,1) / 6.0D+00 s(i+1,4) = ratio * s(i-1,4) + s(i,4) if ( z < 0.5D+00 ) then ratio = s(i,1) * 6.0D+00 * s(i-1,6) * alpha / s(i-1,2) s(i,2) = ratio * s(i-1,3) + s(i,1) + s(i-1,1) s(i,3) = - s(i-1,1) else if ( z == 0.5D+00 ) then ratio = s(i,1) / s(i-1,2) s(i,2) = ratio * s(i-1,3) + s(i,1) + s(i-1,1) s(i,3) = - s(i-1,1) else if ( 0.5D+00 < z ) then factr2 = onemzt * ( alpha * ( onemzt**2 - 1.0D+00 ) + 1.0D+00 ) & / ( alpha * ( onemzt**3 - 1.0D+00 ) + 1.0D+00 ) ratio = factr2 * s(i,1) / s(i-1,2) s(i,2) = ratio * s(i-1,3) + factr2 * s(i-1,1) + s(i,1) s(i,3) = - factr2 * s(i-1,1) end if ! ! At this point, the last two equations read: ! DIAG(I)*XI+ U(I)*XI+1 = R(I) ! -RATIO*DIAG(I)*XI+DIAG(I+1)*XI+1 = R(I+1) ! ! Eliminate XI from the last equation. ! s(i,4) = ratio * s(i-1,4) ratio = - entry / s(i,2) s(i+1,2) = ratio * s(i,3) + s(i+1,2) s(i+1,4) = ratio * s(i,4) + s(i+1,4) ! ! Back substitution. ! s(ntau,4) = s(ntau,4) / s(ntau,2) do while ( 1 < i ) s(i,4) = ( s(i,4) - s(i,3) * s(i+1,4) ) / s(i,2) i = i - 1 end do s(1,4) = ( s(1,4) - s(1,3) * s(2,4) - entry3 * s(3,4) ) / s(1,2) ! ! Construct polynomial pieces. ! break(1) = tau(1) l = 1 do i = 1, ntau - 1 coef(1,l) = gtau(i) coef(3,l) = s(i,4) divdif = ( gtau(i+1) - gtau(i) ) / s(i,1) z = s(i,5) if ( z == 0.0D+00 ) then coef(2,l) = divdif coef(3,l) = 0D+00 coef(4,l) = 0.0D+00 else if ( z < 0.5D+00 ) then zeta = gam * z onemzt = 1.0D+00 - zeta c = s(i+1,4) / 6.0D+00 d = s(i,4) * s(i,6) l = l + 1 del = zeta * s(i,1) break(l) = tau(i) + del zt2 = zeta**2 alpha = alph(onemzt) factor = onemzt**2 * alpha coef(1,l) = gtau(i) + divdif * del & + s(i,1)**2 * ( d * onemzt * ( factor - 1.0D+00 ) & + c * zeta * ( zt2 - 1.0D+00 ) ) coef(2,l) = divdif + s(i,1) * ( d * ( 1.0D+00 - 3.0D+00 * factor ) & + c * ( 3.0D+00 * zt2 - 1.0D+00 ) ) coef(3,l) = 6.0D+00 * ( d * alpha * onemzt + c * zeta ) coef(4,l) = 6.0D+00 * ( c - d * alpha ) / s(i,1) coef(4,l-1) = coef(4,l) & - 6.0D+00 * d * ( 1.0D+00 - alpha ) / ( del * zt2 ) coef(2,l-1) = coef(2,l) - del * ( coef(3,l) & - ( del / 2.0D+00 ) * coef(4,l-1)) else if ( z == 0.5D+00 ) then coef(2,l) = divdif & - s(i,1) * ( 2.0D+00 * s(i,4) + s(i+1,4) ) / 6.0D+00 coef(4,l) = ( s(i+1,4) - s(i,4) ) / s(i,1) else if ( 0.5D+00 <= z ) then onemzt = gam * ( 1.0D+00 - z ) if ( onemzt == 0.0D+00 ) then coef(2,l) = divdif coef(3,l) = 0D+00 coef(4,l) = 0.0D+00 else zeta = 1.0D+00 - onemzt alpha = alph(zeta) c = s(i+1,4) * s(i,6) d = s(i,4) / 6.0D+00 del = zeta * s(i,1) break(l+1) = tau(i) + del coef(2,l) = divdif - s(i,1) * ( 2.0D+00 * d + c ) coef(4,l) = 6.0D+00 * ( c * alpha - d ) / s(i,1) l = l + 1 coef(4,l) = coef(4,l-1) + 6.0D+00 * ( 1.0D+00 - alpha ) * c & / ( s(i,1) * onemzt**3 ) coef(3,l) = coef(3,l-1) + del * coef(4,l-1) coef(2,l) = coef(2,l-1) + del * ( coef(3,l-1) & + ( del / 2.0D+00 ) * coef(4,l-1) ) coef(1,l) = coef(1,l-1) + del * ( coef(2,l-1) & + ( del / 2.0D+00 ) * ( coef(3,l-1) & + ( del / 3.0D+00 ) * coef(4,l-1) ) ) end if end if l = l + 1 break(l) = tau(i+1) end do l = l - 1 k = 4 iflag = 1 return end subroutine titand ( t, g, n ) !*****************************************************************************80 ! !! TITAND represents a temperature-dependent property of titanium. ! ! Discussion: ! ! The data has been used extensively as an example in spline ! approximation with variable knots. ! ! Modified: ! ! 14 February 2007 ! ! Author: ! ! Original FORTRAN77 version by Carl de Boor. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Carl de Boor, ! A Practical Guide to Splines, ! Springer, 2001, ! ISBN: 0387953663, ! LC: QA1.A647.v27. ! ! Parameters: ! ! Output, real ( kind = rk ) T(N), the location of the data points. ! ! Output, real ( kind = rk ) G(N), the value associated with the data points. ! ! Output, integer N, the number of data points, which is 49. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) g(*) integer n real ( kind = rk ) t(*) n = 49 t(1:49) = (/ & 595.0D+00, 605.0D+00, 615.0D+00, 625.0D+00, 635.0D+00, & 645.0D+00, 655.0D+00, 665.0D+00, 675.0D+00, 685.0D+00, & 695.0D+00, 705.0D+00, 715.0D+00, 725.0D+00, 735.0D+00, & 745.0D+00, 755.0D+00, 765.0D+00, 775.0D+00, 785.0D+00, & 795.0D+00, 805.0D+00, 815.0D+00, 825.0D+00, 835.0D+00, & 845.0D+00, 855.0D+00, 865.0D+00, 875.0D+00, 885.0D+00, & 895.0D+00, 905.0D+00, 915.0D+00, 925.0D+00, 935.0D+00, & 945.0D+00, 955.0D+00, 965.0D+00, 975.0D+00, 985.0D+00, & 995.0D+00, 1005.0D+00, 1015.0D+00, 1025.0D+00, 1035.0D+00, & 1045.0D+00, 1055.0D+00, 1065.0D+00, 1075.0D+00 /) g(1:49) = (/ & 0.644D+00, 0.622D+00, 0.638D+00, 0.649D+00, 0.652D+00, & 0.639D+00, 0.646D+00, 0.657D+00, 0.652D+00, 0.655D+00, & 0.644D+00, 0.663D+00, 0.663D+00, 0.668D+00, 0.676D+00, & 0.676D+00, 0.686D+00, 0.679D+00, 0.678D+00, 0.683D+00, & 0.694D+00, 0.699D+00, 0.710D+00, 0.730D+00, 0.763D+00, & 0.812D+00, 0.907D+00, 1.044D+00, 1.336D+00, 1.881D+00, & 2.169D+00, 2.075D+00, 1.598D+00, 1.211D+00, 0.916D+00, & 0.746D+00, 0.672D+00, 0.627D+00, 0.615D+00, 0.607D+00, & 0.606D+00, 0.609D+00, 0.603D+00, 0.601D+00, 0.603D+00, & 0.601D+00, 0.611D+00, 0.601D+00, 0.608D+00 /) return end