src/utils/interpolation.F90

!!
!! Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
!!
!! The University of Notre Dame grants you ("Licensee") a
!! non-exclusive, royalty free, license to use, modify and
!! redistribute this software in source and binary code form, provided
!! that the following conditions are met:
!!
!! 1. Acknowledgement of the program authors must be made in any
!!    publication of scientific results based in part on use of the
!!    program.  An acceptable form of acknowledgement is citation of
!!    the article in which the program was described (Matthew
!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
!!    Parallel Simulation Engine for Molecular Dynamics,"
!!    J. Comput. Chem. 26, pp. 252-271 (2005))
!!
!! 2. Redistributions of source code must retain the above copyright
!!    notice, this list of conditions and the following disclaimer.
!!
!! 3. Redistributions in binary form must reproduce the above copyright
!!    notice, this list of conditions and the following disclaimer in the
!!    documentation and/or other materials provided with the
!!    distribution.
!!
!! This software is provided "AS IS," without a warranty of any
!! kind. All express or implied conditions, representations and
!! warranties, including any implied warranty of merchantability,
!! fitness for a particular purpose or non-infringement, are hereby
!! excluded.  The University of Notre Dame and its licensors shall not
!! be liable for any damages suffered by licensee as a result of
!! using, modifying or distributing the software or its
!! derivatives. In no event will the University of Notre Dame or its
!! licensors be liable for any lost revenue, profit or data, or for
!! direct, indirect, special, consequential, incidental or punitive
!! damages, however caused and regardless of the theory of liability,
!! arising out of the use of or inability to use software, even if the
!! University of Notre Dame has been advised of the possibility of
!! such damages.
!!
!!
!!  interpolation.F90
!!
!!  Created by Charles F. Vardeman II on 03 Apr 2006.
!!
!!  PURPOSE: Generic Spline interplelation routines. These routines assume that we are on a uniform grid for
!!           precomputation of spline parameters.
!!
!! @author Charles F. Vardeman II 
!! @version $Id: interpolation.F90,v 1.3 2006-04-14 21:06:55 chrisfen Exp $


module  INTERPOLATION
  use definitions
  use status
  implicit none
  PRIVATE

  character(len = statusMsgSize) :: errMSG

  type, public :: cubicSpline
     private
     integer :: np = 0
     real(kind=dp) :: dx
     real(kind=dp) :: dx_i
     real (kind=dp), pointer,dimension(:)   :: x => null()
     real (kind=dp), pointer,dimension(:,:) :: c => null()
  end type cubicSpline

  interface newSpline
     module procedure newSpline
  end interface

  public :: deleteSpline

contains


  subroutine newSpline(cs, x, y, yp1, ypn)

    !************************************************************************
    !
    ! newSplineWithoutDerivs solves for slopes defining a cubic spline.
    !
    !  Discussion:
    !
    !    A tridiagonal linear system for the unknown slopes S(I) of
    !    F at x(I), I=1,..., N, is generated and then solved by Gauss
    !    elimination, with S(I) ending up in cs%C(2,I), for all I.
    !
    !  Reference:
    !
    !    Carl DeBoor,
    !    A Practical Guide to Splines,
    !    Springer Verlag.
    !
    !  Parameters:
    !
    !    Input, real x(N), the abscissas or X values of
    !    the data points.  The entries of x are assumed to be
    !    strictly increasing.
    !
    !    Input, real y(I), contains the function value at x(I) for 
    !      I = 1, N.
    !
    !    yp1 contains the slope at x(1) and ypn contains
    !    the slope at x(N).
    !
    !    On output, the intermediate slopes at x(I) have been
    !    stored in cs%C(2,I), for I = 2 to N-1.

    implicit none

    type (cubicSpline), intent(inout) :: cs
    real( kind = DP ), intent(in) :: x(:), y(:)
    real( kind = DP ), intent(in) :: yp1, ypn
    real( kind = DP ) :: g, divdif1, divdif3, dx
    integer :: i, alloc_error, np

    alloc_error = 0

    if (cs%np .ne. 0) then
       call handleWarning("interpolation::newSplineWithoutDerivs", &
            "Type was already created")
       call deleteSpline(cs)
    end if

    ! make sure the sizes match

    if (size(x) .ne. size(y)) then
       call handleError("interpolation::newSplineWithoutDerivs", &
            "Array size mismatch")
    end if

    np = size(x)
    cs%np = np

    allocate(cs%x(np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSplineWithoutDerivs", &
            "Error in allocating storage for x")
    endif

    allocate(cs%c(4,np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSplineWithoutDerivs", &
            "Error in allocating storage for c")
    endif
       
    do i = 1, np
       cs%x(i) = x(i)
       cs%c(1,i) = y(i)       
    enddo

    ! Set the first derivative of the function to the second coefficient of 
    ! each of the endpoints

    cs%c(2,1) = yp1
    cs%c(2,np) = ypn
    

    !
    !  Set up the right hand side of the linear system.
    !
    do i = 2, cs%np - 1
       cs%c(2,i) = 3.0_DP * ( &
            (x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
            (x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
    end do
    !
    !  Set the diagonal coefficients.
    !
    cs%c(4,1) = 1.0_DP
    do i = 2, cs%np - 1 
       cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
    end do
    cs%c(4,cs%np) = 1.0_DP
    !
    !  Set the off-diagonal coefficients.
    !
    cs%c(3,1) = 0.0_DP
    do i = 2, cs%np
       cs%c(3,i) = x(i) - x(i-1)
    end do
    !
    !  Forward elimination.
    !
    do i = 2, cs%np - 1
       g = -cs%c(3,i+1) / cs%c(4,i-1)
       cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
       cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
    end do
    !
    !  Back substitution for the interior slopes.
    !
    do i = cs%np - 1, 2, -1
       cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
    end do
    !
    !  Now compute the quadratic and cubic coefficients used in the 
    !  piecewise polynomial representation.
    !
    do i = 1, cs%np - 1
       dx = x(i+1) - x(i)
       divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
       divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
       cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
       cs%c(4,i) = divdif3 / ( dx * dx )
    end do

    cs%c(3,cs%np) = 0.0_DP
    cs%c(4,cs%np) = 0.0_DP

    cs%dx = dx
    cs%dx_i = 1.0_DP / dx
    return
  end subroutine newSplineWithoutDerivs

  subroutine deleteSpline(this)

    type(cubicSpline) :: this
    
    if(associated(this%x)) then
       deallocate(this%x)
       this%x => null()
    end if
    if(associated(this%c)) then
       deallocate(this%c)
       this%c => null()
    end if
    
    this%np = 0
    
  end subroutine deleteSpline

  subroutine lookup_nonuniform_spline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    newSpline must be called first, to set up the
    !    spline data from the raw function and derivative data.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Reference:
    !
    !    Conte and de Boor,
    !    Algorithm PCUBIC,
    !    Elementary Numerical Analysis, 
    !    1973, page 234.
    !
    !  Parameters:
    !
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) :: dx
    integer :: i, j
    !
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.
    !
    j = cs%np - 1

    do i = 1, cs%np - 2

       if ( xval < cs%x(i+1) ) then
          j = i
          exit
       end if

    end do
    !
    !  Evaluate the cubic polynomial.
    !
    dx = xval - cs%x(j)

    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
    
    return
  end subroutine lookup_nonuniform_spline

  subroutine lookup_uniform_spline(cs, xval, yval)
    
    !*************************************************************************
    !
    ! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
    !
    !  Discussion:
    !
    !    newSpline must be called first, to set up the
    !    spline data from the raw function and derivative data.
    !
    !  Modified:
    !
    !    06 April 1999
    !
    !  Reference:
    !
    !    Conte and de Boor,
    !    Algorithm PCUBIC,
    !    Elementary Numerical Analysis, 
    !    1973, page 234.
    !
    !  Parameters:
    !
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) :: dx
    integer :: i, j
    !
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.

    j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))

    dx = xval - cs%x(j)

    yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
    
    return
  end subroutine lookup_uniform_spline
  
end module INTERPOLATION
Revision:	2710
Committed:	Fri Apr 14 21:06:55 2006 UTC (18 years, 2 months ago) by chrisfen
File size:	9519 byte(s)
Log Message:	Wiped the spline with derivatives and corrected a boundary comparision error by eliminating the check. This algorithm REQUIRES known first derivatives at the endpoints.
#	User	Rev	Content
1	gezelter	2708	!!
2			!! Copyright (c) 2006 The University of Notre Dame. All Rights Reserved.
3			!!
4			!! The University of Notre Dame grants you ("Licensee") a
5			!! non-exclusive, royalty free, license to use, modify and
6			!! redistribute this software in source and binary code form, provided
7			!! that the following conditions are met:
8			!!
9			!! 1. Acknowledgement of the program authors must be made in any
10			!! publication of scientific results based in part on use of the
11			!! program. An acceptable form of acknowledgement is citation of
12			!! the article in which the program was described (Matthew
13			!! A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			!! J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			!! Parallel Simulation Engine for Molecular Dynamics,"
16			!! J. Comput. Chem. 26, pp. 252-271 (2005))
17			!!
18			!! 2. Redistributions of source code must retain the above copyright
19			!! notice, this list of conditions and the following disclaimer.
20			!!
21			!! 3. Redistributions in binary form must reproduce the above copyright
22			!! notice, this list of conditions and the following disclaimer in the
23			!! documentation and/or other materials provided with the
24			!! distribution.
25			!!
26			!! This software is provided "AS IS," without a warranty of any
27			!! kind. All express or implied conditions, representations and
28			!! warranties, including any implied warranty of merchantability,
29			!! fitness for a particular purpose or non-infringement, are hereby
30			!! excluded. The University of Notre Dame and its licensors shall not
31			!! be liable for any damages suffered by licensee as a result of
32			!! using, modifying or distributing the software or its
33			!! derivatives. In no event will the University of Notre Dame or its
34			!! licensors be liable for any lost revenue, profit or data, or for
35			!! direct, indirect, special, consequential, incidental or punitive
36			!! damages, however caused and regardless of the theory of liability,
37			!! arising out of the use of or inability to use software, even if the
38			!! University of Notre Dame has been advised of the possibility of
39			!! such damages.
40			!!
41			!!
42			!! interpolation.F90
43			!!
44			!! Created by Charles F. Vardeman II on 03 Apr 2006.
45			!!
46			!! PURPOSE: Generic Spline interplelation routines. These routines assume that we are on a uniform grid for
47			!! precomputation of spline parameters.
48			!!
49			!! @author Charles F. Vardeman II
50	chrisfen	2710	!! @version $Id: interpolation.F90,v 1.3 2006-04-14 21:06:55 chrisfen Exp $
51	gezelter	2708
52
53			module INTERPOLATION
54			use definitions
55			use status
56			implicit none
57			PRIVATE
58
59			character(len = statusMsgSize) :: errMSG
60
61			type, public :: cubicSpline
62			private
63			integer :: np = 0
64			real(kind=dp) :: dx
65			real(kind=dp) :: dx_i
66			real (kind=dp), pointer,dimension(:) :: x => null()
67	gezelter	2709	real (kind=dp), pointer,dimension(:,:) :: c => null()
68	gezelter	2708	end type cubicSpline
69
70			interface newSpline
71	chrisfen	2710	module procedure newSpline
72	gezelter	2708	end interface
73
74			public :: deleteSpline
75
76			contains
77
78
79	chrisfen	2710	subroutine newSpline(cs, x, y, yp1, ypn)
80	gezelter	2708
81			!************************************************************************
82			!
83			! newSplineWithoutDerivs solves for slopes defining a cubic spline.
84			!
85			! Discussion:
86			!
87			! A tridiagonal linear system for the unknown slopes S(I) of
88			! F at x(I), I=1,..., N, is generated and then solved by Gauss
89			! elimination, with S(I) ending up in cs%C(2,I), for all I.
90			!
91			! Reference:
92			!
93			! Carl DeBoor,
94			! A Practical Guide to Splines,
95			! Springer Verlag.
96			!
97			! Parameters:
98			!
99			! Input, real x(N), the abscissas or X values of
100	chrisfen	2710	! the data points. The entries of x are assumed to be
101	gezelter	2708	! strictly increasing.
102			!
103			! Input, real y(I), contains the function value at x(I) for
104			! I = 1, N.
105			!
106			! yp1 contains the slope at x(1) and ypn contains
107			! the slope at x(N).
108			!
109			! On output, the intermediate slopes at x(I) have been
110			! stored in cs%C(2,I), for I = 2 to N-1.
111
112			implicit none
113
114			type (cubicSpline), intent(inout) :: cs
115			real( kind = DP ), intent(in) :: x(:), y(:)
116			real( kind = DP ), intent(in) :: yp1, ypn
117			real( kind = DP ) :: g, divdif1, divdif3, dx
118			integer :: i, alloc_error, np
119
120			alloc_error = 0
121
122			if (cs%np .ne. 0) then
123			call handleWarning("interpolation::newSplineWithoutDerivs", &
124			"Type was already created")
125			call deleteSpline(cs)
126			end if
127
128			! make sure the sizes match
129
130			if (size(x) .ne. size(y)) then
131			call handleError("interpolation::newSplineWithoutDerivs", &
132			"Array size mismatch")
133			end if
134
135			np = size(x)
136			cs%np = np
137
138			allocate(cs%x(np), stat=alloc_error)
139			if(alloc_error .ne. 0) then
140			call handleError("interpolation::newSplineWithoutDerivs", &
141			"Error in allocating storage for x")
142			endif
143
144			allocate(cs%c(4,np), stat=alloc_error)
145			if(alloc_error .ne. 0) then
146			call handleError("interpolation::newSplineWithoutDerivs", &
147			"Error in allocating storage for c")
148			endif
149
150			do i = 1, np
151			cs%x(i) = x(i)
152			cs%c(1,i) = y(i)
153			enddo
154
155	chrisfen	2710	! Set the first derivative of the function to the second coefficient of
156			! each of the endpoints
157	gezelter	2708
158	chrisfen	2710	cs%c(2,1) = yp1
159			cs%c(2,np) = ypn
160
161
162	gezelter	2708	!
163			! Set up the right hand side of the linear system.
164			!
165			do i = 2, cs%np - 1
166			cs%c(2,i) = 3.0_DP * ( &
167			(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
168			(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
169			end do
170			!
171			! Set the diagonal coefficients.
172			!
173			cs%c(4,1) = 1.0_DP
174			do i = 2, cs%np - 1
175			cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
176			end do
177	gezelter	2709	cs%c(4,cs%np) = 1.0_DP
178	gezelter	2708	!
179			! Set the off-diagonal coefficients.
180			!
181			cs%c(3,1) = 0.0_DP
182			do i = 2, cs%np
183			cs%c(3,i) = x(i) - x(i-1)
184			end do
185			!
186			! Forward elimination.
187			!
188			do i = 2, cs%np - 1
189			g = -cs%c(3,i+1) / cs%c(4,i-1)
190			cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
191			cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
192			end do
193			!
194			! Back substitution for the interior slopes.
195			!
196			do i = cs%np - 1, 2, -1
197			cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
198			end do
199			!
200			! Now compute the quadratic and cubic coefficients used in the
201			! piecewise polynomial representation.
202			!
203			do i = 1, cs%np - 1
204			dx = x(i+1) - x(i)
205			divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
206			divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
207			cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
208			cs%c(4,i) = divdif3 / ( dx * dx )
209			end do
210
211	gezelter	2709	cs%c(3,cs%np) = 0.0_DP
212			cs%c(4,cs%np) = 0.0_DP
213	gezelter	2708
214			cs%dx = dx
215	gezelter	2709	cs%dx_i = 1.0_DP / dx
216	gezelter	2708	return
217			end subroutine newSplineWithoutDerivs
218
219			subroutine deleteSpline(this)
220
221			type(cubicSpline) :: this
222
223			if(associated(this%x)) then
224			deallocate(this%x)
225			this%x => null()
226			end if
227			if(associated(this%c)) then
228			deallocate(this%c)
229			this%c => null()
230			end if
231
232			this%np = 0
233
234			end subroutine deleteSpline
235
236			subroutine lookup_nonuniform_spline(cs, xval, yval)
237
238			!*************************************************************************
239			!
240			! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
241			!
242			! Discussion:
243			!
244			! newSpline must be called first, to set up the
245			! spline data from the raw function and derivative data.
246			!
247			! Modified:
248			!
249			! 06 April 1999
250			!
251			! Reference:
252			!
253			! Conte and de Boor,
254			! Algorithm PCUBIC,
255			! Elementary Numerical Analysis,
256			! 1973, page 234.
257			!
258			! Parameters:
259			!
260			implicit none
261
262			type (cubicSpline), intent(in) :: cs
263			real( kind = DP ), intent(in) :: xval
264			real( kind = DP ), intent(out) :: yval
265	gezelter	2709	real( kind = DP ) :: dx
266	gezelter	2708	integer :: i, j
267			!
268			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
269			! or is nearest to xval.
270			!
271			j = cs%np - 1
272
273			do i = 1, cs%np - 2
274
275			if ( xval < cs%x(i+1) ) then
276			j = i
277			exit
278			end if
279
280			end do
281			!
282			! Evaluate the cubic polynomial.
283			!
284			dx = xval - cs%x(j)
285
286			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
287
288			return
289			end subroutine lookup_nonuniform_spline
290
291			subroutine lookup_uniform_spline(cs, xval, yval)
292
293			!*************************************************************************
294			!
295			! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
296			!
297			! Discussion:
298			!
299			! newSpline must be called first, to set up the
300			! spline data from the raw function and derivative data.
301			!
302			! Modified:
303			!
304			! 06 April 1999
305			!
306			! Reference:
307			!
308			! Conte and de Boor,
309			! Algorithm PCUBIC,
310			! Elementary Numerical Analysis,
311			! 1973, page 234.
312			!
313			! Parameters:
314			!
315			implicit none
316
317			type (cubicSpline), intent(in) :: cs
318			real( kind = DP ), intent(in) :: xval
319			real( kind = DP ), intent(out) :: yval
320	gezelter	2709	real( kind = DP ) :: dx
321	gezelter	2708	integer :: i, j
322			!
323			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
324			! or is nearest to xval.
325
326	gezelter	2709	j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
327	gezelter	2708
328			dx = xval - cs%x(j)
329
330			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
331
332			return
333			end subroutine lookup_uniform_spline
334
335			end module INTERPOLATION