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 interpolation 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.4 2006-04-14 21:49:54 gezelter Exp $


module  INTERPOLATION
  use definitions
  use status
  implicit none
  PRIVATE

  character(len = statusMsgSize) :: errMSG

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

  public :: newSpline
  public :: deleteSpline
  public :: lookup_spline
  public :: lookup_uniform_spline
  public :: lookup_nonuniform_spline
  
contains
  

  subroutine newSpline(cs, x, y, yp1, ypn, isUniform)
    
    !************************************************************************
    !
    ! newSpline 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.
    !
    !    Input, real yp1 contains the slope at x(1) 
    !    Input, real ypn contains the slope at x(N)
    !
    !    On output, the slopes at x(I) have been stored in 
    !               cs%C(2,I), for I = 1 to N.

    implicit none

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

    alloc_error = 0

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

    ! make sure the sizes match

    np = size(x)

    if ( size(y) .ne. np ) then
       call handleError("interpolation::newSpline", &
            "Array size mismatch")
    end if
    
    cs%np = np
    cs%isUniform = isUniform

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

    allocate(cs%c(4,np), stat=alloc_error)
    if(alloc_error .ne. 0) then
       call handleError("interpolation::newSpline", &
            "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_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

  subroutine lookup_spline(cs, xval, yval)

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(inout) :: xval
    real( kind = DP ), intent(inout) :: yval
    
    if (cs%isUniform) then
       call lookup_uniform_spline(cs, xval, yval)
    else
       call lookup_nonuniform_spline(cs, xval, yval)
    endif

    return
  end subroutine lookup_spline
  
end module INTERPOLATION
Revision:	2711
Committed:	Fri Apr 14 21:49:54 2006 UTC (18 years, 2 months ago) by gezelter
File size:	10005 byte(s)
Log Message:	added some logical flags
#	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	gezelter	2711	!! PURPOSE: Generic Spline interpolation routines. These routines
47			!! assume that we are on a uniform grid for precomputation of
48			!! spline parameters.
49	gezelter	2708	!!
50			!! @author Charles F. Vardeman II
51	gezelter	2711	!! @version $Id: interpolation.F90,v 1.4 2006-04-14 21:49:54 gezelter Exp $
52	gezelter	2708
53
54			module INTERPOLATION
55			use definitions
56			use status
57			implicit none
58			PRIVATE
59
60			character(len = statusMsgSize) :: errMSG
61
62			type, public :: cubicSpline
63			private
64	gezelter	2711	logical :: isUniform = .false.
65	gezelter	2708	integer :: np = 0
66			real(kind=dp) :: dx_i
67			real (kind=dp), pointer,dimension(:) :: x => null()
68	gezelter	2709	real (kind=dp), pointer,dimension(:,:) :: c => null()
69	gezelter	2708	end type cubicSpline
70
71	gezelter	2711	public :: newSpline
72	gezelter	2708	public :: deleteSpline
73	gezelter	2711	public :: lookup_spline
74			public :: lookup_uniform_spline
75			public :: lookup_nonuniform_spline
76
77	gezelter	2708	contains
78	gezelter	2711
79	gezelter	2708
80	gezelter	2711	subroutine newSpline(cs, x, y, yp1, ypn, isUniform)
81
82	gezelter	2708	!************************************************************************
83			!
84	gezelter	2711	! newSpline solves for slopes defining a cubic spline.
85	gezelter	2708	!
86			! Discussion:
87			!
88			! A tridiagonal linear system for the unknown slopes S(I) of
89			! F at x(I), I=1,..., N, is generated and then solved by Gauss
90			! elimination, with S(I) ending up in cs%C(2,I), for all I.
91			!
92			! Reference:
93			!
94			! Carl DeBoor,
95			! A Practical Guide to Splines,
96			! Springer Verlag.
97			!
98			! Parameters:
99			!
100			! Input, real x(N), the abscissas or X values of
101	chrisfen	2710	! the data points. The entries of x are assumed to be
102	gezelter	2708	! strictly increasing.
103			!
104			! Input, real y(I), contains the function value at x(I) for
105			! I = 1, N.
106			!
107	gezelter	2711	! Input, real yp1 contains the slope at x(1)
108			! Input, real ypn contains the slope at x(N)
109	gezelter	2708	!
110	gezelter	2711	! On output, the slopes at x(I) have been stored in
111			! cs%C(2,I), for I = 1 to N.
112	gezelter	2708
113			implicit none
114
115			type (cubicSpline), intent(inout) :: cs
116			real( kind = DP ), intent(in) :: x(:), y(:)
117			real( kind = DP ), intent(in) :: yp1, ypn
118	gezelter	2711	logical, intent(in) :: isUniform
119	gezelter	2708	real( kind = DP ) :: g, divdif1, divdif3, dx
120			integer :: i, alloc_error, np
121
122			alloc_error = 0
123
124			if (cs%np .ne. 0) then
125	gezelter	2711	call handleWarning("interpolation::newSpline", &
126			"cubicSpline struct was already created")
127	gezelter	2708	call deleteSpline(cs)
128			end if
129
130			! make sure the sizes match
131
132	gezelter	2711	np = size(x)
133
134			if ( size(y) .ne. np ) then
135			call handleError("interpolation::newSpline", &
136	gezelter	2708	"Array size mismatch")
137			end if
138	gezelter	2711
139	gezelter	2708	cs%np = np
140	gezelter	2711	cs%isUniform = isUniform
141	gezelter	2708
142			allocate(cs%x(np), stat=alloc_error)
143			if(alloc_error .ne. 0) then
144	gezelter	2711	call handleError("interpolation::newSpline", &
145	gezelter	2708	"Error in allocating storage for x")
146			endif
147
148			allocate(cs%c(4,np), stat=alloc_error)
149			if(alloc_error .ne. 0) then
150	gezelter	2711	call handleError("interpolation::newSpline", &
151	gezelter	2708	"Error in allocating storage for c")
152			endif
153
154			do i = 1, np
155			cs%x(i) = x(i)
156			cs%c(1,i) = y(i)
157			enddo
158
159	chrisfen	2710	! Set the first derivative of the function to the second coefficient of
160			! each of the endpoints
161	gezelter	2708
162	chrisfen	2710	cs%c(2,1) = yp1
163			cs%c(2,np) = ypn
164
165	gezelter	2708	!
166			! Set up the right hand side of the linear system.
167			!
168	gezelter	2711
169	gezelter	2708	do i = 2, cs%np - 1
170			cs%c(2,i) = 3.0_DP * ( &
171			(x(i) - x(i-1)) * (cs%c(1,i+1) - cs%c(1,i)) / (x(i+1) - x(i)) + &
172			(x(i+1) - x(i)) * (cs%c(1,i) - cs%c(1,i-1)) / (x(i) - x(i-1)))
173			end do
174	gezelter	2711
175	gezelter	2708	!
176			! Set the diagonal coefficients.
177			!
178			cs%c(4,1) = 1.0_DP
179			do i = 2, cs%np - 1
180			cs%c(4,i) = 2.0_DP * ( x(i+1) - x(i-1) )
181			end do
182	gezelter	2709	cs%c(4,cs%np) = 1.0_DP
183	gezelter	2708	!
184			! Set the off-diagonal coefficients.
185			!
186			cs%c(3,1) = 0.0_DP
187			do i = 2, cs%np
188			cs%c(3,i) = x(i) - x(i-1)
189			end do
190			!
191			! Forward elimination.
192			!
193			do i = 2, cs%np - 1
194			g = -cs%c(3,i+1) / cs%c(4,i-1)
195			cs%c(4,i) = cs%c(4,i) + g * cs%c(3,i-1)
196			cs%c(2,i) = cs%c(2,i) + g * cs%c(2,i-1)
197			end do
198			!
199			! Back substitution for the interior slopes.
200			!
201			do i = cs%np - 1, 2, -1
202			cs%c(2,i) = ( cs%c(2,i) - cs%c(3,i) * cs%c(2,i+1) ) / cs%c(4,i)
203			end do
204			!
205			! Now compute the quadratic and cubic coefficients used in the
206			! piecewise polynomial representation.
207			!
208			do i = 1, cs%np - 1
209			dx = x(i+1) - x(i)
210			divdif1 = ( cs%c(1,i+1) - cs%c(1,i) ) / dx
211			divdif3 = cs%c(2,i) + cs%c(2,i+1) - 2.0_DP * divdif1
212			cs%c(3,i) = ( divdif1 - cs%c(2,i) - divdif3 ) / dx
213			cs%c(4,i) = divdif3 / ( dx * dx )
214			end do
215
216	gezelter	2709	cs%c(3,cs%np) = 0.0_DP
217			cs%c(4,cs%np) = 0.0_DP
218	gezelter	2708
219	gezelter	2709	cs%dx_i = 1.0_DP / dx
220	gezelter	2711
221	gezelter	2708	return
222			end subroutine newSplineWithoutDerivs
223
224			subroutine deleteSpline(this)
225
226			type(cubicSpline) :: this
227
228			if(associated(this%x)) then
229			deallocate(this%x)
230			this%x => null()
231			end if
232			if(associated(this%c)) then
233			deallocate(this%c)
234			this%c => null()
235			end if
236
237			this%np = 0
238
239			end subroutine deleteSpline
240
241			subroutine lookup_nonuniform_spline(cs, xval, yval)
242
243			!*************************************************************************
244			!
245			! lookup_nonuniform_spline evaluates a piecewise cubic Hermite interpolant.
246			!
247			! Discussion:
248			!
249			! newSpline must be called first, to set up the
250			! spline data from the raw function and derivative data.
251			!
252			! Modified:
253			!
254			! 06 April 1999
255			!
256			! Reference:
257			!
258			! Conte and de Boor,
259			! Algorithm PCUBIC,
260			! Elementary Numerical Analysis,
261			! 1973, page 234.
262			!
263			! Parameters:
264			!
265			implicit none
266
267			type (cubicSpline), intent(in) :: cs
268			real( kind = DP ), intent(in) :: xval
269			real( kind = DP ), intent(out) :: yval
270	gezelter	2709	real( kind = DP ) :: dx
271	gezelter	2708	integer :: i, j
272			!
273			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
274			! or is nearest to xval.
275			!
276			j = cs%np - 1
277
278			do i = 1, cs%np - 2
279
280			if ( xval < cs%x(i+1) ) then
281			j = i
282			exit
283			end if
284
285			end do
286			!
287			! Evaluate the cubic polynomial.
288			!
289			dx = xval - cs%x(j)
290
291			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
292
293			return
294			end subroutine lookup_nonuniform_spline
295
296			subroutine lookup_uniform_spline(cs, xval, yval)
297
298			!*************************************************************************
299			!
300			! lookup_uniform_spline evaluates a piecewise cubic Hermite interpolant.
301			!
302			! Discussion:
303			!
304			! newSpline must be called first, to set up the
305			! spline data from the raw function and derivative data.
306			!
307			! Modified:
308			!
309			! 06 April 1999
310			!
311			! Reference:
312			!
313			! Conte and de Boor,
314			! Algorithm PCUBIC,
315			! Elementary Numerical Analysis,
316			! 1973, page 234.
317			!
318			! Parameters:
319			!
320			implicit none
321
322			type (cubicSpline), intent(in) :: cs
323			real( kind = DP ), intent(in) :: xval
324			real( kind = DP ), intent(out) :: yval
325	gezelter	2709	real( kind = DP ) :: dx
326	gezelter	2708	integer :: i, j
327			!
328			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
329			! or is nearest to xval.
330
331	gezelter	2709	j = MAX(1, MIN(cs%np, idint((xval-cs%x(1)) * cs%dx_i) + 1))
332	gezelter	2708
333			dx = xval - cs%x(j)
334
335			yval = cs%c(1,j) + dx * ( cs%c(2,j) + dx * ( cs%c(3,j) + dx * cs%c(4,j) ) )
336
337			return
338			end subroutine lookup_uniform_spline
339	gezelter	2711
340			subroutine lookup_spline(cs, xval, yval)
341
342			type (cubicSpline), intent(in) :: cs
343			real( kind = DP ), intent(inout) :: xval
344			real( kind = DP ), intent(inout) :: yval
345
346			if (cs%isUniform) then
347			call lookup_uniform_spline(cs, xval, yval)
348			else
349			call lookup_nonuniform_spline(cs, xval, yval)
350			endif
351
352			return
353			end subroutine lookup_spline
354	gezelter	2708
355			end module INTERPOLATION