DISPATCH
lagrange_mod.f90
1 !===============================================================================
2 !> Module for Lagrange interpolation
3 !===============================================================================
5  USE io_mod
6  USE trace_mod
7  implicit none
8  private
9  type, public:: lagrange_t
10  integer:: verbose=0
11  integer:: order=2
12  logical:: initialized=.false.
13  contains
14  procedure interpolation4
15  procedure interpolation8
16  generic:: interpolation => interpolation4, interpolation8
17  procedure sequence4
18  procedure sequence8
19  generic:: sequence => sequence4, sequence8
20  procedure interval4
21  procedure interval8
22  generic:: interval => interval4, interval8
23  procedure weights4
24  procedure weights8
25  generic:: weights => weights4, weights8
26  procedure deriv_weights4
27  procedure deriv_weights8
28  generic:: deriv_weights => deriv_weights4, deriv_weights8
29  procedure deriv2_weights4
30  procedure deriv2_weights8
31  generic:: deriv2_weights => deriv2_weights4, deriv2_weights8
32  procedure sequence_weights4
33  procedure sequence_weights8
34  generic:: sequence_weights => sequence_weights4, sequence_weights8
35  procedure deriv_sequence_weights4
36  procedure deriv_sequence_weights8
37  generic:: deriv_sequence_weights => deriv_sequence_weights4, deriv_sequence_weights8
38  procedure deriv2_sequence_weights4
39  procedure deriv2_sequence_weights8
40  generic:: deriv2_sequence_weights => deriv2_sequence_weights4, deriv2_sequence_weights8
41  procedure interpolate3d
42  procedure trilinear3d
43  procedure test
44  procedure debug
45  procedure init
46  procedure interpolate1d4
47  procedure interpolate1d8
48  generic:: interpolate1d => interpolate1d4, interpolate1d8
49  end type
50  type(lagrange_t), public:: lagrange
51 CONTAINS
52 
53 !> -----------------------------------------------------------------------------
54 !> Classical Lagrange interpolation
55 !> -----------------------------------------------------------------------------
56 FUNCTION interpolation4 (self, x0, x, y) RESULT(f)
57  class(lagrange_t):: self
58  real(kind=4), dimension(:):: x, y
59  real(kind=4):: x0, f, g
60  integer:: n, i, j
61  !.............................................................................
62  f = 0.0
63  n = size(x)
64  do i=1,n
65  g = y(i)
66  do j=1,n
67  g = merge(g,g*(x0-x(j))/(x(i)-x(j)),j==i)
68  end do
69  f = f + g
70  end do
71 END FUNCTION interpolation4
72 FUNCTION interpolation8 (self, x0, x, y) RESULT(f)
73  class(lagrange_t):: self
74  real(kind=8), dimension(:):: x, y
75  real(kind=8):: x0, f, g
76  integer:: n, i, j
77  !.............................................................................
78  f = 0.0
79  n = size(x)
80  do i=1,n
81  g = y(i)
82  do j=1,n
83  g = merge(g,g*(x0-x(j))/(x(i)-x(j)),j==i)
84  end do
85  f = f + g
86  end do
87 END FUNCTION interpolation8
88 
89 !> -----------------------------------------------------------------------------
90 !> Classical Lagrange interpolation weights. w(i) is the contribution factor
91 !> of a value y(i) to the function value at x0, given that there are mesh points
92 !> at x = [x(j), j=1,n], so the function value f(x0) = sum(y(i)*w(i))
93 !> -----------------------------------------------------------------------------
94 FUNCTION weights4 (self, x0, x) RESULT(w)
95  class(lagrange_t):: self
96  real(kind=4), dimension(:):: x
97  real(kind=4), dimension(size(x)):: w
98  real(kind=4):: x0, g
99  integer:: n, i, j
100  !.............................................................................
101  n = size(x)
102  do i=1,n
103  w(i) = 1.0
104  do j=1,n
105  if (j==i) cycle
106  if (x(j)==x(i)) then
107  w = 0.0
108  w(1) = 1.0
109  return
110  end if
111  w(i) = w(i)*(x0-x(j))/(x(i)-x(j))
112  end do
113  end do
114 END FUNCTION weights4
115 FUNCTION weights8 (self, x0, x) RESULT(w)
116  class(lagrange_t):: self
117  real(kind=8), dimension(:):: x
118  real(kind=8), dimension(size(x)):: w
119  real(kind=8):: x0, g
120  integer:: n, i, j
121  !.............................................................................
122  n = size(x)
123  do i=1,n
124  w(i) = 1.0
125  do j=1,n
126  if (j==i) cycle
127  if (x(j)==x(i)) then
128  w = 0d0
129  w(1) = 1d0
130  return
131  end if
132  w(i) = w(i)*(x0-x(j))/(x(i)-x(j))
133  end do
134  end do
135 END FUNCTION weights8
136 
137 !> -----------------------------------------------------------------------------
138 !> Classical Lagrange interpolation weights. w(i) is the contribution factor
139 !> of a value y(i) to the derivative value at x0, given that there are mesh
140 !> points at x = [x(j), j=1,n], so the function value f(x0) = sum(y(i)*w(i))
141 !> -----------------------------------------------------------------------------
142 FUNCTION deriv_weights4 (self, x0, x) RESULT(w)
143  class(lagrange_t):: self
144  real(kind=4), dimension(:):: x
145  real(kind=4), dimension(size(x)):: w
146  real(kind=4):: x0, nomin, denom, g
147  integer:: n, i, j, k
148  !.............................................................................
149  n = size(x)
150  do i=1,n
151  denom = 1.0_4
152  nomin = 0.0_4
153  do j=1,n
154  if (j==i) cycle
155  if (x(i)==x(j)) then
156  w = 0.0
157  return
158  end if
159  g = 1.
160  do k=1,n
161  if (k==j .or. k==i) cycle
162  g = g*(x0-x(k))
163  end do
164  nomin = nomin + g
165  denom = denom*(x(i)-x(j))
166  end do
167  w(i) = nomin/denom
168  end do
169 END FUNCTION deriv_weights4
170 FUNCTION deriv_weights8 (self, x0, x) RESULT(w)
171  class(lagrange_t):: self
172  real(kind=8), dimension(:):: x
173  real(kind=8), dimension(size(x)):: w
174  real(kind=8):: x0, nomin, denom, g
175  integer:: n, i, j, k
176  !.............................................................................
177  n = size(x)
178  do i=1,n
179  denom = 1.0_8
180  nomin = 0.0_8
181  do j=1,n
182  if (j==i) cycle
183  if (x(i)==x(j)) then
184  w = 0.0_8
185  return
186  end if
187  g = 1.0_8
188  do k=1,n
189  if (k==j .or. k==i) cycle
190  g = g*(x0-x(k))
191  end do
192  nomin = nomin + g
193  denom = denom*(x(i)-x(j))
194  end do
195  w(i) = nomin/denom
196  end do
197 END FUNCTION deriv_weights8
198 
199 !> -----------------------------------------------------------------------------
200 !> Classical Lagrange interpolation weights. w(i) is the contribution factor
201 !> of a value y(i) to the 2nd derivative value at x0, given that there are mesh
202 !> points at x = [x(j), j=1,n], so the function value f(x0) = sum(y(i)*w(i))
203 !> -----------------------------------------------------------------------------
204 FUNCTION deriv2_weights4 (self, x0, x) RESULT(w)
205  class(lagrange_t):: self
206  real(kind=4), dimension(:):: x
207  real(kind=4), dimension(size(x)):: w
208  real(kind=4):: x0, nomin, denom, g
209  integer:: n, i, j, k, l
210  !.............................................................................
211  n = size(x)
212  do i=1,n
213  denom = 1.0_4
214  nomin = 0.0_4
215  do j=1,n
216  if (j==i) cycle
217  if (x(i)==x(j)) then
218  w = 0.0
219  return
220  end if
221  do k=1,n
222  if (k==j .or. k==i) cycle
223  g = 1.
224  do l=1,n
225  if (l==k .or. l==j .or. l==i) cycle
226  g = g*(x0-x(l))
227  end do
228  nomin = nomin + g
229  end do
230  denom = denom*(x(i)-x(j))
231  end do
232  w(i) = nomin/denom
233  end do
234 END FUNCTION deriv2_weights4
235 FUNCTION deriv2_weights8 (self, x0, x) RESULT(w)
236  class(lagrange_t):: self
237  real(kind=8), dimension(:):: x
238  real(kind=8), dimension(size(x)):: w
239  real(kind=8):: x0, nomin, denom, g
240  integer:: n, i, j, k, l
241  !.............................................................................
242  n = size(x)
243  do i=1,n
244  denom = 1.0_8
245  nomin = 0.0_8
246  do j=1,n
247  if (j==i) cycle
248  if (x(i)==x(j)) then
249  w = 0.0_8
250  return
251  end if
252  do k=1,n
253  if (k == j .or. k == i) cycle
254  g = 1.0_8
255  do l=1,n
256  if (l == k .or. l == j .or. l == i) cycle
257  g = g*(x0-x(l))
258  end do
259  nomin = nomin + g
260  end do
261  denom = denom*(x(i)-x(j))
262  end do
263  w(i) = nomin/denom
264  end do
265 END FUNCTION deriv2_weights8
266 
267 !> -----------------------------------------------------------------------------
268 !> Prepare to use Lagrange interpolation of given order to interpolate to
269 !> coordinate x0 in a sequence that may be longer than needed
270 !> -----------------------------------------------------------------------------
271 SUBROUTINE interval4 (self, x0, x, i0, i2, o, order)
272  class(lagrange_t):: self
273  integer, optional:: order
274  real(kind=4), dimension(:):: x
275  real(kind=4), dimension(size(x)):: y
276  real(kind=4):: x0, f
277  integer:: i0, i1, i2, o, n
278  !.............................................................................
279  ! Default order = 3, or n-1 if that is smaller
280  n = size(x)
281  o = 3
282  if (present(order)) o = order
283  o = min(o,n-1)
284  ! Search for identical initial points
285  i1 = 1
286  do while (x(i1)==x(1) .and. i1<n)
287  i1 = i1+1
288  end do
289  ! All points are equal to x(1)
290  if (x(i1)==x(1)) then
291  i0 = 1
292  i2 = 1
293  o = 0
294  return
295  ! Only the last point is not equal to x(1)
296  else if (i1==n) then
297  i0 = n-1
298  i2 = n
299  o = 1
300  return
301  end if
302  ! Normal search
303  do while (x(i1)<x0 .and. i1<n)
304  i1 = i1+1
305  end do
306  if (i1 < 3 .or. i1 > n-1) then
307  o = min(o,2)
308  end if
309 
310  if (i1 < 2+o/2) then
311  ! We are at the start of the interval, and need 1 point before i1
312  i0 = i1 - (o+1)/2
313  if (i0 < 1) then
314  i0 = 1
315  end if
316  i2 = i0+o
317  if (i2 > n) then
318  i2 = n
319  o = min(o,i2-i0)
320  end if
321  else
322  ! We are in the middle or end of the interval
323  i2 = i1 + (o-1)/2
324  ! Shift left if not enough points
325  if (i2 > n) then
326  i2 = n
327  end if
328  i0 = i2 - o
329  if (i0 < 1) then
330  i0 = 1
331  i2 = min(i2,n)
332  o = min(o,i2-i0)
333  end if
334  end if
335  if (self%verbose>0) then
336  print*,'lagrange_interval: i0, i2, o =', i0, i2, o
337  end if
338 END SUBROUTINE interval4
339 SUBROUTINE interval8 (self, x0, x, i0, i2, o, order)
340  class(lagrange_t):: self
341  integer, optional:: order
342  real(kind=8), dimension(:):: x
343  real(kind=8), dimension(size(x)):: y
344  real(kind=8):: x0, f
345  integer:: i0, i1, i2, o, n
346  !.............................................................................
347  ! Default order = 3, or n-1 if that is smaller
348  n = size(x)
349  o = 3
350  if (present(order)) o = order
351  o = min(o,n-1)
352  ! Search for identical initial points
353  i1 = 1
354  do while (x(i1)==x(1) .and. i1<n)
355  i1 = i1+1
356  end do
357  ! All points are equal to x(1)
358  if (x(i1)==x(1)) then
359  i0 = 1
360  i2 = 1
361  o = 0
362  return
363  ! Only the last point is not equal to x(1)
364  else if (i1==n) then
365  i0 = n-1
366  i2 = n
367  o = 1
368  return
369  end if
370  ! Normal search
371  do while (x(i1)<x0 .and. i1<n)
372  i1 = i1+1
373  end do
374  if (i1 < 3 .or. i1 > n-1) then
375  o = min(o,2)
376  end if
377 
378  if (i1 < 2+o/2) then
379  ! We are at the start of the interval, and need 1 point before i1
380  i0 = i1 - (o+1)/2
381  if (i0 < 1) then
382  i0 = 1
383  end if
384  i2 = i0+o
385  if (i2 > n) then
386  i2 = n
387  o = min(o,i2-i0)
388  end if
389  else
390  ! We are in the middle or end of the interval
391  i2 = i1 + (o-1)/2
392  ! Shift left if not enough points
393  if (i2 > n) then
394  i2 = n
395  end if
396  i0 = i2 - o
397  if (i0 < 1) then
398  i0 = 1
399  i2 = min(i2,n)
400  o = min(o,i2-i0)
401  end if
402  end if
403  if (self%debug(1)) then
404  print '("lagrange_interval: i0,i2:",2i3,3x,"o:",i2,3x,"x0:",f10.6,3x,"x:",10f10.6)', &
405  i0,i2, o, x0, x
406  end if
407 END SUBROUTINE interval8
408 
409 !> -----------------------------------------------------------------------------
410 !> Use Lagrange interpolation of given order to interpolate to coordinate x0
411 !> in a sequence that may be longer than needed
412 !> -----------------------------------------------------------------------------
413 FUNCTION sequence4 (self, x0, x, y, order, correct) RESULT(f)
414  class(lagrange_t):: self
415  integer, optional:: order
416  real(kind=4), optional:: correct
417  real(kind=4), dimension(:):: x, y
418  real(kind=4):: x0, f
419  integer:: i0, i2, o, n
420  !.............................................................................
421  n = size(x)
422  call self%interval (x0, x, i0, i2, o, order)
423  if (self%verbose>0) then
424  print'(" i0,i2,order=",3i3,4x,"x0,x,y =",1p,18g12.4)', &
425  i0, i2, o, x0, x(i0:i2), y(i0:i2)
426  end if
427  f = self%interpolation (x0, x(i0:i2), y(i0:i2))
428  if (self%verbose>0) then
429  if (present(correct)) then
430  print'(" i0,i2,order=",3i3,4x,"x0,x,y,f,correct=",1p,18g12.4)', &
431  i0, i2, o, x0, x(i0:i2), y(i0:i2), f, correct
432  else
433  print'(" i0,i2,order=",3i3,4x,"x0,x,y,f=",1p,18g12.4)', &
434  i0, i2, o, x0, x(i0:i2), y(i0:i2), f
435  end if
436  end if
437 END FUNCTION sequence4
438 FUNCTION sequence8 (self, x0, x, y, order, correct) RESULT(f)
439  class(lagrange_t):: self
440  integer, optional:: order
441  real(kind=8), optional:: correct
442  real(kind=8), dimension(:):: x, y
443  real(kind=8):: x0, f
444  integer:: i0, i2, o, n
445  !.............................................................................
446  n = size(x)
447  call self%interval (x0, x, i0, i2, o, order)
448  if (self%verbose>0) then
449  print'(" i0,i2,order=",3i3,4x,"x0,x,y =",1p,18g12.4)', &
450  i0, i2, o, x0, x(i0:i2), y(i0:i2)
451  end if
452  f = self%interpolation (x0, x(i0:i2), y(i0:i2))
453  if (self%verbose>0) then
454  if (present(correct)) then
455  print'(" i0,i2,order=",3i3,4x,"x0,x,y,f,correct=",1p,18g12.4)', &
456  i0, i2, o, x0, x(i0:i2), y(i0:i2), f, correct
457  else
458  print'(" i0,i2,order=",3i3,4x,"x0,x,y,f=",1p,18g12.4)', &
459  i0, i2, o, x0, x(i0:i2), y(i0:i2), f
460  end if
461  end if
462 END FUNCTION sequence8
463 
464 !> -----------------------------------------------------------------------------
465 !> Use Lagrange interpolation of given order to get weights for interpolate to
466 !> coordinate x0 in a sequence of increasing x that may be longer than needed
467 !> -----------------------------------------------------------------------------
468 SUBROUTINE sequence_weights4 (self, x0, x, i0, i2, w, order)
469  class(lagrange_t):: self
470  integer, optional:: order
471  real(kind=4), dimension(:):: x, w
472  real(kind=4):: x0
473  integer:: i0, i2, o, n
474  !.............................................................................
475  n = size(x)
476  call self%interval (x0, x, i0, i2, o, order)
477  w = 0d0
478  if (o==0) then
479  w(1) = 1.0
480  else
481  w(1:i2-i0+1) = self%weights (x0, x(i0:i2))
482  end if
483 END SUBROUTINE sequence_weights4
484 SUBROUTINE sequence_weights8 (self, x0, x, i0, i2, w, order)
485  class(lagrange_t):: self
486  integer, optional:: order
487  real(kind=8), dimension(:):: x, w
488  real(kind=8):: x0
489  integer:: i0, i2, o, n
490  !.............................................................................
491  n = size(x)
492  call self%interval (x0, x, i0, i2, o, order)
493  w = 0d0
494  if (o==0) then
495  w(1) = 1.0
496  else
497  w(1:i2-i0+1) = self%weights (x0, x(i0:i2))
498  end if
499 END SUBROUTINE sequence_weights8
500 
501 !> -----------------------------------------------------------------------------
502 !> Use Lagrange interpolation of given order to get weights for interpolate to
503 !> coordinate x0 in a sequence of increasing x that may be longer than needed
504 !> -----------------------------------------------------------------------------
505 SUBROUTINE deriv_sequence_weights4 (self, x0, x, i0, i2, w, order)
506  class(lagrange_t):: self
507  integer, optional:: order
508  real(kind=4), dimension(:):: x, w
509  real(kind=4):: x0
510  integer:: i0, i2, o, n
511  !.............................................................................
512  n = size(x)
513  call self%interval (x0, x, i0, i2, o, order)
514  w = 0d0
515  w(1:i2-i0+1) = self%deriv_weights (x0, x(i0:i2))
516 END SUBROUTINE deriv_sequence_weights4
517 SUBROUTINE deriv_sequence_weights8 (self, x0, x, i0, i2, w, order)
518  class(lagrange_t):: self
519  integer, optional:: order
520  real(kind=8), dimension(:):: x, w
521  real(kind=8):: x0
522  integer:: i0, i2, o, n
523  !.............................................................................
524  n = size(x)
525  call self%interval (x0, x, i0, i2, o, order)
526  w = 0d0
527  w(1:i2-i0+1) = self%deriv_weights (x0, x(i0:i2))
528 END SUBROUTINE deriv_sequence_weights8
529 
530 !> -----------------------------------------------------------------------------
531 !> Use Lagrange interpolation of given order to get weights for interpolate to
532 !> coordinate x0 in a sequence of increasing x that may be longer than needed
533 !> -----------------------------------------------------------------------------
534 SUBROUTINE deriv2_sequence_weights4 (self, x0, x, i0, i2, w, order)
535  class(lagrange_t):: self
536  integer, optional:: order
537  real(kind=4), dimension(:):: x, w
538  real(kind=4):: x0
539  integer:: i0, i2, o, n
540  !.............................................................................
541  n = size(x)
542  call self%interval (x0, x, i0, i2, o, order)
543  w = 0d0
544  w(1:i2-i0+1) = self%deriv2_weights (x0, x(i0:i2))
545 END SUBROUTINE deriv2_sequence_weights4
546 SUBROUTINE deriv2_sequence_weights8 (self, x0, x, i0, i2, w, order)
547  class(lagrange_t):: self
548  integer, optional:: order
549  real(kind=8), dimension(:):: x, w
550  real(kind=8):: x0
551  integer:: i0, i2, o, n
552  !.............................................................................
553  n = size(x)
554  call self%interval (x0, x, i0, i2, o, order)
555  w = 0d0
556  w(1:i2-i0+1) = self%deriv2_weights (x0, x(i0:i2))
557 END SUBROUTINE deriv2_sequence_weights8
558 
559 !> -----------------------------------------------------------------------------
560 !> Interpolate in 3D from fi to fo
561 !> -----------------------------------------------------------------------------
562 SUBROUTINE interpolate3d (self, xi, yi, zi, fi, xo, yo, zo, fo, order)
563  class(lagrange_t):: self
564  real, dimension(:):: xi, yi, zi, xo, yo, zo
565  real, dimension(:,:,:):: fi, fo
566  integer:: order
567  !.............................................................................
568  real, dimension(:,:,:), pointer:: f1, f2
569  integer:: nix, niy, niz, nox, noy, noz, ix, iy, iz
570  !-----------------------------------------------------------------------------
571  nix = size(xi)
572  niy = size(yi)
573  niz = size(zi)
574  nox = size(xo)
575  noy = size(yo)
576  noz = size(zo)
577  allocate (f1(nix,niy,noz))
578  do iz=1,noz
579  do iy=1,niy
580  do ix=1,nix
581  f1(ix,iy,iz) = self%sequence (zo(iz), zi, fi(ix,iy,:), order)
582  end do
583  end do
584  end do
585  allocate (f2(nix,noy,noz))
586  do iz=1,noz
587  do iy=1,noy
588  do ix=1,nix
589  f2(ix,iy,iz) = self%sequence (yo(iy), yi, f1(ix,:,iz), order)
590  end do
591  end do
592  end do
593  deallocate (f1)
594  do iz=1,noz
595  do iy=1,noy
596  do ix=1,nox
597  fo(ix,iy,iz) = self%sequence (xo(ix), xi, f2(:,iy,iz), order)
598  end do
599  end do
600  end do
601  deallocate (f2)
602 END SUBROUTINE interpolate3d
603 
604 !> -----------------------------------------------------------------------------
605 !> Interpolate in 3D from fi to fo
606 !> -----------------------------------------------------------------------------
607 SUBROUTINE trilinear3d (self, xi, yi, zi, fi, xo, yo, zo, fo, order)
608  class(lagrange_t):: self
609  real, dimension(:):: xi, yi, zi, xo, yo, zo
610  real, dimension(:,:,:):: fi, fo
611  integer:: order
612  !.............................................................................
613  integer:: nix, niy, niz, nox, noy, noz, ix, iy, iz, jx ,jy, jz
614  real:: px, py, pz, qx, qy, qz, dx, dy, dz
615  !-----------------------------------------------------------------------------
616  nix = size(xi)
617  niy = size(yi)
618  niz = size(zi)
619  nox = size(xo)
620  noy = size(yo)
621  noz = size(zo)
622  dx = (xi(nix)-xi(1))/(nix-1)
623  dy = (yi(niy)-yi(1))/(niy-1)
624  dz = (zi(niz)-zi(1))/(niz-1)
625  do iz=1,noz
626  do jz=1,niz-2
627  if (zi(jz+1)>zo(iz)) exit
628  end do
629  pz = (zo(iz)-zi(jz))/(zi(jz+1)-zi(jz))
630  pz = max(0.0,min(1.0,pz))
631  qz = 1.0-pz
632  do iy=1,noy
633  py = (yo(iy)-yi(1))/dy + 1.
634  jy = py
635  jy = max(1,min(niy-1,jy))
636  py = max(0.0,min(1.0,py -jy))
637  qy = 1.0-py
638  do ix=1,nox
639  px = (xo(ix)-xi(1))/dx + 1.
640  jx = px
641  jx = max(1,min(nix-1,jx))
642  px = max(0.0,min(1.0,px -jx))
643  qx = 1.0-px
644  fo(ix,iy,iz) = &
645  qz*(qy*(qx*fi(jx ,jy , jz ) + px*fi(jx+1,jy , jz )) + &
646  py*(qx*fi(jx ,jy+1, jz ) + px*fi(jx+1,jy+1, jz ))) + &
647  pz*(qy*(qx*fi(jx ,jy , jz+1) + px*fi(jx+1,jy , jz+1)) + &
648  py*(qx*fi(jx ,jy+1, jz+1) + px*fi(jx+1,jy+1, jz+1)))
649  end do
650  end do
651  end do
652 END SUBROUTINE trilinear3d
653 
654 !> -----------------------------------------------------------------------------
655 !> Test for a polynomial of order 2
656 !> -----------------------------------------------------------------------------
657 SUBROUTINE test (self)
658  class(lagrange_t):: self
659  integer, parameter:: n=5
660  integer:: order, i, i0, i1
661  real(8), dimension(n):: x, y, w
662  real(8):: x0, y0
663  !-----------------------------------------------------------------------------
664  if (.not. io_unit%master) return
665  print *,'--------------- lagrange_t%test ------------------'
666  x = 0.0
667  y = 1.0
668  order = 1
669  print '(5x,a,11x,a,11x,a,4x,a)', 'x', 'y', 'y==1', 'order'
670  do i=1,n+1
671  x0 = i-0.5
672  call self%sequence_weights (x0, x, i0, i1, w, order)
673  y0 = sum(y(i0:i1)*w(1:1+i1-i0))
674  print '(3g12.4,i6)', x0, y0, x0**order, i1-i0
675  end do
676  y = 0.0
677  x(n) = 1.0
678  y(n) = 1.0
679  print '(5x,a,11x,a,11x,a,4x,a)', 'x', 'y', 'x-4', 'order'
680  do i=1,n+1
681  x0 = i-0.5
682  call self%sequence_weights (x0, x, i0, i1, w, order)
683  y0 = sum(y(i0:i1)*w(1:1+i1-i0))
684  print '(3g12.4,i6)', x0, y0, x0**order, i1-i0
685  end do
686  order = 2
687  do i=1,n
688  x(i) = i
689  y(i) = x(i)**order
690  end do
691  print '(5x,a,11x,a,11x,a,4x,a)', 'x', 'y', 'x**2', 'order'
692  do i=1,n+1
693  x0 = i-0.5
694  call self%sequence_weights (x0, x, i0, i1, w, order)
695  y0 = sum(y(i0:i1)*w(1:1+i1-i0))
696  print '(3g12.4,i6)', x0, y0, x0**order, i1-i0
697  end do
698  print *,'--------------------------------------------------'
699 END SUBROUTINE test
700 
701 !> -----------------------------------------------------------------------------
702 !> -----------------------------------------------------------------------------
703 FUNCTION debug (self, verbose)
704  class(lagrange_t):: self
705  integer:: verbose
706  logical:: debug
707  !-----------------------------------------------------------------------------
708  if (.not.self%initialized) call self%init
709  debug = (self%verbose >= verbose)
710 END FUNCTION debug
711 
712 !> -----------------------------------------------------------------------------
713 !> Initialize (only) this instance of lagrange
714 !> -----------------------------------------------------------------------------
715 SUBROUTINE init (self)
716  class(lagrange_t):: self
717  integer:: iostat
718  integer:: verbose=0, order=2
719  namelist /lagrange_params/ verbose, order
720  !-----------------------------------------------------------------------------
721  call trace%begin ('lagrange_t%init')
722  !$omp critical (input_cr)
723  if (.not.self%initialized) then
724  self%initialized = .true.
725  rewind(io_unit%input)
726  read (io_unit%input, lagrange_params, iostat=iostat)
727  if (io_unit%master) write (io_unit%output, lagrange_params)
728  self%verbose = verbose
729  self%order = order
730  end if
731  !$omp end critical (input_cr)
732  call trace%end ()
733 END SUBROUTINE init
734 
735 !> =============================================================================
736 !> Use default 1-D Lagrange interpolation
737 !> =============================================================================
738 FUNCTION interpolate1d4 (self,x0, n, x, y) RESULT(f)
739  class(lagrange_t):: self
740  integer:: n
741  real(kind=4):: x0, x(n), y(n), f
742  !.............................................................................
743  if (x0 < x(1)) then
744  f = y(1)
745  else if (x0 > x(n)) then
746  f = y(n)
747  else
748  f = self%sequence (x0, x, y)
749  end if
750 END FUNCTION interpolate1d4
751 
752 FUNCTION interpolate1d8 (self,x0, n, x, y) RESULT(f)
753  class(lagrange_t):: self
754  integer:: n
755  real(kind=8):: x0, x(n), y(n), f
756  !.............................................................................
757  if (x0 < x(1)) then
758  f = y(1)
759  else if (x0 > x(n)) then
760  f = y(n)
761  else
762  f = self%sequence (x0, x, y)
763  end if
764 END FUNCTION interpolate1d8
765 
766 END MODULE
Module for Lagrange interpolation.
Definition: lagrange_mod.f90:4
Definition: io_mod.f90:4