!> !> @file splopt.f90 !> !> @brief !> !> @copyright !> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) !> SPC (Swiss Plasma Center) !> !> SPClibs is free software: you can redistribute it and/or modify it under !> the terms of the GNU Lesser General Public License as published by the Free !> Software Foundation, either version 3 of the License, or (at your option) !> any later version. !> !> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY !> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS !> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. !> !> You should have received a copy of the GNU Lesser General Public License !> along with this program. If not, see . !> !> @author !> (in alphabetical order) !> @author Trach-Minh Tran !> subroutine splopt ( tau, n, k, scrtch, t, iflag ) !************************************************************************* ! !! 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), ..., T(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) + (-)**(n-k-j)*x(j) , j=1,...,n-k, ! with X the solution of the linear system ! C * X = D. ! ! Here, c(i,j)=b(i)(xi(j)), all j, with b(i) the i-th b-spline of ! order K for the knot sequence TAU, all i, and D is the vector ! given by d(i)=sum( -a(j) , j=i,...,n )*(TAU(i+k)-TAU(i))/k, all i, ! with a(i)=sum ( (-)**(n-k-j)*b(i,k+1,tau)(xi(j)) , j=1,...,n-k ) ! for i=1,...,n-1, and a(n)=-.5 . ! ! See chapter XIII of text and references there for a derivation. ! ! The first guess for t(k+j) is (TAU(j+1)+...+TAU(j+k-1))/(k-1) . ! iteration terminates if max(abs(x(j))) < t o l , with ! TOL = t o l r t e *(TAU(n)-TAU(1))/(n-k) , ! or else after NEWTMX iterations , currently, ! newtmx, tolrte / 10, .000001 ! ! Reference: ! ! Carl DeBoor, ! A Practical Guide to Splines, ! Springer Verlag. ! ! Parameters: ! ! Input, real ( kind = 8 ) TAU(N), the interpolation points. ! assumed to be nondecreasing, with tau(i) < tau(i+k),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 = 8 ) SCRTCH((N-K)*(2*K+3)+5*K+3). The various ! contents are specified in the text below . ! ! Output, real ( kind = 8 ) T(N+K), the optimal knots ready for ! use in optimal recovery. specifically, t(1)=... = t(k) = ! tau(1) and t(n+1)=... = t(n+k) = tau(n) , while the n-k ! interior knots t(k+1), ..., t(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 k integer n real ( kind = 8 ) del real ( kind = 8 ) delmax real ( kind = 8 ) floatk integer i integer id integer iflag integer index integer j integer kp1 integer kpkm1 integer kpn integer l 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 = 8 ) scrtch((n-k)*(2*k+3)+5*k+3) real ( kind = 8 ) t(n+k) real ( kind = 8 ) tau(n) real ( kind = 8 ) sign real ( kind = 8 ) signst real ( kind = 8 ) sum1 real ( kind = 8 ) tol real ( kind = 8 ), parameter :: tolrte = 0.000001D+00 real ( kind = 8 ) xij nmk = n - k if ( n < k ) then write ( *, * ) ' ' write ( *, * ) 'SPLOPT - Fatal error!' write ( *, * ) ' N < K, N = ',n,' K = ',k iflag = 2 return end if if ( n == k ) then do i = 1, k t(i) = tau(1) t(n+i) = tau(n) end do return end if if ( k <= 2 ) then write(*,*)' ' write(*,*)'SPLOPT - Fatal error!' write(*,*)' K < 2, K = ',k iflag = 2 stop 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. ! do j = 1, k scrtch(j) = tau(1) scrtch(kpn+j) = tau(n) end do do j = 1, n scrtch(k+j) = tau(j) end do ! ! First guess for scrtch (.+nx) = xi . ! scrtch(nx) = tau(1) scrtch(nmk+1+nx) = tau(n) do j = 1, nmk sum1 = 0.0D+00 do l = 1, k-1 sum1 = sum1 + tau(j+l) end do scrtch(j+nx) = sum1 / real ( k - 1, kind = 8 ) 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 = 8 ) ! ! Start the Newton step. ! compute the 2k-1 bands of the matrix c and store in scrtch(.+nc), ! and compute the vector scrtch(.+na)=-a. ! 100 continue do i = 1, lenw scrtch(i+nc) = 0.0D+00 end do do i = 2, n scrtch(i-1+na) = 0.0D+00 end do sign = signst left = kp1 do j = 1, nmk xij = scrtch(j+nx) 130 continue if ( xij < scrtch(left+1) ) then go to 140 end if left = left+1 if ( left < kpn ) then go to 130 end if left = left-1 140 continue 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 2k-1 bands of c in the 2k-1 r o w s of ! the work array w, and w in turn is stored in s c r t c h , ! with w(1,1)=scrtch(1+nc). ! ! also, c being of order n-k, we would want ! 1 <= left-2k+ll .le. n-k or ! llmin=2k-left <= ll .le. 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 ( *, * ) ' ' write ( *, * ) 'SPLOPT - Fatal error!' write ( *, * ) ' Matrix C is not invertible.' stop 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 go to 230 end if del = max ( del, ( scrtch(i-1+nx) - scrtch(i+nx) ) / 3.0D+00 ) go to 240 230 del = min (del,(scrtch(i+1+nx)-scrtch(i+nx))/3.0D+00 ) 240 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 go to 270 end if newton = newton + 1 if ( newton <= newtmx ) then go to 100 end if write ( *, * ) ' ' write ( *, * ) 'SPLOPT - Warning!' write ( *, * ) ' No convergence. Number of Newton steps was ', newtmx 270 continue do i = 1, nmk t(k+i) = scrtch(i+nx) end do 290 continue do i=1,k t(i)=tau(1) t(n+i)=tau(n) end do return ! 310 iflag=2 ! ! return end