lag_interp.F90

Go to the documentation of this file.

!> Fortran module to prepare for Lagrange polynomial interpolation.
!> based on GSI: lagmod.f90
 
module lag_interp_mod

use kinds, only: kind_real
use gnssro_mod_constants,only: zero, one

implicit none

! set default to private
private
! set subroutines to public
public :: lag_interp_const
public :: lag_interp_const_tl
public :: lag_interp_const_ad
public :: lag_interp_weights
public :: lag_interp_weights_tl
public :: lag_interp_weights_ad
public :: lag_interp_smthweights
public :: lag_interp_smthweights_tl
public :: lag_interp_smthweights_ad


contains

subroutine lag_interp_const(q,x,n)
! Precompute the N constant denominator factors of the N-point 
! Lagrange polynomial interpolation formula.
!
! input argument list:
!   X -    The N abscissae.
!   N -    The number of points involved.
!
! output argument list:
!   Q -    The N denominator constants.
IMPLICIT NONE

INTEGER          ,INTENT(in   ) :: n
REAL(kind_real),INTENT(  out) :: q(n)
REAL(kind_real),INTENT(in   ) :: x(n)
!-----------------------------------------------------------------------------
INTEGER                      :: i,j
!=============================================================================
DO i=1,n
   q(i)=one
   DO j=1,n
      IF(j /= i)q(i)=q(i)/(x(i)-x(j))
   ENDDO
ENDDO
end subroutine lag_interp_const


!============================================================================
subroutine lag_interp_const_tl(q,q_TL,x,x_TL,n)
!============================================================================
IMPLICIT NONE

INTEGER          ,INTENT(in   ) :: n !number of points involved
REAL(kind_real),DIMENSION(n),INTENT(  out) :: q,q_tl
REAL(kind_real),DIMENSION(n),INTENT(in   ) :: x,x_tl
!-----------------------------------------------------------------------------
INTEGER                      :: i,j
REAL(kind_real)                         :: rat
!=============================================================================
DO i=1,n
   q(i)=one
   q_tl(i)=zero
   DO j=1,n
      IF(j /= i) THEN
         rat=one/(x(i)-x(j))
         q_tl(i)=(q_tl(i)-q(i)*(x_tl(i)-x_tl(j))*rat)*rat
         q(i)=q(i)*rat
      ENDIF
   ENDDO
ENDDO
end subroutine lag_interp_const_tl


!============================================================================
subroutine lag_interp_const_ad(q_AD,x,x_AD,n)
!============================================================================
IMPLICIT NONE

INTEGER          ,INTENT(in   ) :: n
REAL(kind_real),DIMENSION(n),INTENT(in   ) :: x
REAL(kind_real),DIMENSION(n),INTENT(inout) :: x_ad
REAL(kind_real),DIMENSION(n),INTENT(inout) :: q_ad
!-----------------------------------------------------------------------------
INTEGER                      :: i,j
REAL(kind_real),DIMENSION(n)            :: q
REAL(kind_real),DIMENSION(n,n)          :: jac
!=============================================================================
call lag_interp_const(q,x,n)
jac=zero
DO i=1,n
   DO j=1,n
      IF(j /= i) THEN
         jac(j,i)=q(i)/(x(i)-x(j))
         jac(i,i)=jac(i,i)-jac(j,i)
      ENDIF
   ENDDO
ENDDO
x_ad=x_ad+matmul(jac,q_ad)
end subroutine lag_interp_const_ad

!============================================================================
subroutine lag_interp_weights(x,xt,q,w,dw,n)
! Construct the Lagrange weights and their derivatives when 
! target abscissa is known and denominators Q have already 
! been precomputed
!
! input argument list:
!   X   - Grid abscissae
!   XT  - Target abscissa
!   Q   - Q factors (denominators of the Lagrange weight formula)
!   N   - Number of grid points involved in the interpolation
!
! output argument list:
!   W   - Lagrange weights
!   DW  - Derivatives, dW/dX, of Lagrange weights W
!============================================================================
IMPLICIT NONE

INTEGER          ,INTENT(IN   ) :: n
REAL(kind_real)             ,INTENT(IN   ) :: xt
REAL(kind_real),INTENT(IN   ) :: x(n),q(n)
REAL(kind_real),INTENT(  OUT) :: w(n),dw(n)
!-----------------------------------------------------------------------------
REAL(kind_real)            :: d(n),pa(n),pb(n),dpa(n),dpb(n)
INTEGER                      :: j
!============================================================================
pa(1)=one
dpa(1)=zero
do j=1,n-1
   d(j)=xt-x(j)
   pa(j+1)=pa(j)*d(j)
   dpa(j+1)=dpa(j)*d(j)+pa(j)
enddo
d(n)=xt-x(n)

pb(n)=one
dpb(n)=zero
do j=n,2,-1
   pb(j-1)=pb(j)*d(j)
   dpb(j-1)=dpb(j)*d(j)+pb(j)
enddo
do j=1,n
   w(j)= pa(j)*pb(j)*q(j)
   dw(j)=(pa(j)*dpb(j)+dpa(j)*pb(j))*q(j)
enddo
end subroutine lag_interp_weights


!============================================================================
subroutine lag_interp_weights_tl(x,x_TL,xt,q,q_TL,w,w_TL,dw,dw_TL,n)
!============================================================================
IMPLICIT NONE

INTEGER          ,INTENT(IN   ) :: n
REAL(kind_real)             ,INTENT(IN   ) :: xt
REAL(kind_real),DIMENSION(n),INTENT(IN   ) :: x,q,x_tl,q_tl
REAL(kind_real),DIMENSION(n),INTENT(  OUT) :: w,dw,w_tl,dw_tl
!-----------------------------------------------------------------------------
REAL(kind_real),DIMENSION(n)            :: d,pa,pb,dpa,dpb
REAL(kind_real),DIMENSION(n)            :: d_tl,pa_tl,pb_tl,dpa_tl,dpb_tl
INTEGER                      :: j
!============================================================================
pa(1)=one
dpa(1)=zero
pa_tl(1)=zero
dpa_tl(1)=zero

do j=1,n-1
   d(j)=xt-x(j)
   d_tl(j)=-x_tl(j)
   pa(j+1)=pa(j)*d(j)
   pa_tl(j+1)=pa_tl(j)*d(j)+pa(j)*d_tl(j)
   dpa(j+1)=dpa(j)*d(j)+pa(j)
   dpa_tl(j+1)=dpa_tl(j)*d(j)+dpa(j)*d_tl(j)+pa_tl(j)
enddo
d(n)=xt-x(n)
d_tl(n)=-x_tl(n)

pb(n)=one
dpb(n)=zero
pb_tl(n)=zero
dpb_tl(n)=zero
do j=n,2,-1
   pb(j-1)=pb(j)*d(j)
   pb_tl(j-1)=pb_tl(j)*d(j)+pb(j)*d_tl(j)
   dpb(j-1)=dpb(j)*d(j)+pb(j)
   dpb_tl(j-1)=dpb_tl(j)*d(j)+dpb(j)*d_tl(j)+pb_tl(j)
enddo
do j=1,n
   w(j)= pa(j)*pb(j)*q(j)
   dw(j)=(pa(j)*dpb(j)+dpa(j)*pb(j))*q(j)
   w_tl(j)=(pa_tl(j)*pb(j)+pa(j)*pb_tl(j))*q(j)+pa(j)*pb(j)*q_tl(j)
   dw_tl(j)=(pa_tl(j)*dpb(j)+pa(j)*dpb_tl(j)+dpa_tl(j)*pb(j)+dpa(j)*pb_tl(j))*q(j)+ &
            (pa(j)*dpb(j)+dpa(j)*pb(j))*q_tl(j)
enddo
end subroutine lag_interp_weights_tl


!============================================================================
subroutine lag_interp_weights_ad(x,x_AD,xt,q,q_AD,w_AD,dw_AD,n)

!============================================================================
IMPLICIT NONE

INTEGER          ,INTENT(IN   ) :: n
REAL(kind_real)             ,INTENT(IN   ) :: xt
REAL(kind_real),DIMENSION(n),INTENT(IN   ) :: x,q
REAL(kind_real),DIMENSION(n),INTENT(INOUT) :: q_ad,x_ad,w_ad,dw_ad
!-----------------------------------------------------------------------------
REAL(kind_real),DIMENSION(n)            :: d,pa,pb,dpa,dpb
REAL(kind_real),DIMENSION(n)            :: d_ad,pa_ad,pb_ad,dpa_ad,dpb_ad
INTEGER                      :: j
!============================================================================

pa_ad=zero; dpb_ad=zero; dpa_ad=zero; pb_ad=zero
d_ad=zero

! passive variables
pa(1)=one
dpa(1)=zero
do j=1,n-1
   d(j)=xt-x(j)
   pa(j+1)=pa(j)*d(j)
   dpa(j+1)=dpa(j)*d(j)+pa(j)
enddo
d(n)=xt-x(n)
pb(n)=one
dpb(n)=zero
do j=n,2,-1
   pb(j-1)=pb(j)*d(j)
   dpb(j-1)=dpb(j)*d(j)+pb(j)
enddo
!
do j=n,1,-1
   pa_ad(j)=pa_ad(j)+ dpb(j)*q(j)*dw_ad(j)
   dpb_ad(j)=dpb_ad(j)+ pa(j)*q(j)*dw_ad(j)
   dpa_ad(j)=dpa_ad(j)+ pb(j)*q(j)*dw_ad(j)
   pb_ad(j)=pb_ad(j)+ dpa(j)*q(j)*dw_ad(j)
   q_ad(j)=q_ad(j)+(pa(j)*dpb(j)+dpa(j)*pb(j))*dw_ad(j)
   pa_ad(j)=pa_ad(j)+ pb(j)*q(j)*w_ad(j)
   pb_ad(j)=pb_ad(j)+ pa(j)*q(j)*w_ad(j)
   q_ad(j)=q_ad(j)+ pa(j)*pb(j)*w_ad(j)
end do
do j=2,n
   dpb_ad(j)=dpb_ad(j)+d(j)*dpb_ad(j-1)
   d_ad(j)=d_ad(j)+dpb(j)*dpb_ad(j-1)
   pb_ad(j)=pb_ad(j)       +dpb_ad(j-1)
   pb_ad(j)=pb_ad(j)+d(j)*pb_ad(j-1)
   d_ad(j)=d_ad(j)+pb(j)*pb_ad(j-1)
enddo

dpb_ad(n)=zero ! not sure about it
pb_ad(n) =zero ! not sure about it

x_ad(n)=x_ad(n)-d_ad(n)
do j=n-1,1,-1
   dpa_ad(j)=dpa_ad(j)+d(j)*dpa_ad(j+1)
   d_ad(j)=d_ad(j)+dpa(j)*dpa_ad(j+1)
   pa_ad(j)=pa_ad(j)       +dpa_ad(j+1)
   pa_ad(j)=pa_ad(j)+d(j)*pa_ad(j+1)
   d_ad(j)=d_ad(j)+pa(j)*pa_ad(j+1)
   x_ad(j)=x_ad(j)       -d_ad(j  )
enddo

end subroutine lag_interp_weights_ad


!============================================================================
subroutine lag_interp_smthweights(x,xt,aq,bq,w,dw,n)
! Construct weights and their derivatives for interpolation 
! to a given target based on a linearly weighted mixture of 
! the pair of Lagrange interpolators which omit the respective 
! end points of the source nodes provided. The number of source 
! points provided must be even and the nominal target interval 
! is the unique central one. The linear weighting pair is 
! determined by the relative location of the target within 
! this central interval, or else the extreme values, 0 and 1, 
! when target lies outside this interval.  The objective is to 
! provide an interpolator whose derivative is continuous.
!============================================================================
IMPLICIT NONE

INTEGER            ,INTENT(IN   ) :: n
REAL(kind_real)               ,INTENT(IN   ) :: xt
REAL(kind_real),INTENT(IN   ) :: x(n)
REAL(kind_real),INTENT(IN   ) :: aq(n-1),bq(n-1)
REAL(kind_real),INTENT(  OUT) :: w(n),dw(n)
!-----------------------------------------------------------------------------
REAL(kind_real)               :: aw(n),bw(n),daw(n),dbw(n)
REAL(kind_real)               :: xa,xb,dwb,wb
INTEGER                       :: na
!============================================================================
CALL lag_interp_weights(x(1:n-1),xt,aq,aw(1:n-1),daw(1:n-1),n-1)
CALL lag_interp_weights(x(2:n  ),xt,bq,bw(2:n  ),dbw(2:n  ),n-1)
aw(n)=zero
daw(n)=zero
bw(1)=zero
dbw(1)=zero
na=n/2
IF(na*2 /= n)stop 'In lag_interp_smthWeights; n must be even'
xa =x(na     )
xb =x(na+1)
dwb=one/(xb-xa)
wb =(xt-xa)*dwb
IF    (wb>one )THEN
   wb =one
   dwb=zero
ELSEIF(wb<zero)THEN
   wb =zero
   dwb=zero
ENDIF
bw =bw -aw
dbw=dbw-daw

w =aw +wb*bw
dw=daw+wb*dbw+dwb*bw
end subroutine lag_interp_smthweights


!============================================================================
subroutine lag_interp_smthweights_tl(x,x_TL,xt,aq,aq_TL,bq,bq_TL,dw,dw_TL,n)
!============================================================================
IMPLICIT NONE

INTEGER            ,INTENT(IN   ) :: n
REAL(kind_real)               ,INTENT(IN   ) :: xt
REAL(kind_real),DIMENSION(n)  ,INTENT(IN   ) :: x,x_tl
REAL(kind_real),DIMENSION(n-1),INTENT(IN   ) :: aq,bq,aq_tl,bq_tl
REAL(kind_real),DIMENSION(n)  ,INTENT(  OUT) :: dw_tl,dw
!-----------------------------------------------------------------------------
REAL(kind_real),DIMENSION(n)               :: aw,bw,daw,dbw
REAL(kind_real),DIMENSION(n)               :: aw_tl,bw_tl,daw_tl,dbw_tl
REAL(kind_real)                            :: xa,xb,dwb,wb
REAL(kind_real)                            :: xa_tl,xb_tl,dwb_tl,wb_tl
INTEGER                         :: na
!============================================================================
CALL lag_interp_weights_tl(x(1:n-1),x_tl(1:n-1),xt,aq,aq_tl,aw(1:n-1),aw_tl(1:n-1),&
             daw(1:n-1),daw_tl(1:n-1),n-1)
CALL lag_interp_weights_tl(x(2:n     ),x_tl(2:n     ),xt,bq,bq_tl,bw(2:n     ),bw_tl(2:n     ),&
             dbw(2:n     ),dbw_tl(2:n     ),n-1)
aw(n) =zero
daw(n)=zero
bw(1) =zero
dbw(1)=zero
!
aw_tl(n) =zero
daw_tl(n)=zero
bw_tl(1) =zero
dbw_tl(1)=zero
na=n/2
IF(na*2 /= n)stop 'In lag_interp_smthWeights; n must be even'
xa   =x(na)
xa_tl=x_tl(na)
xb   =x(na+1)
xb_tl=x_tl(na+1)
dwb   = one          /(xb-xa)
dwb_tl=-(xb_tl-xa_tl)/(xb-xa)**2
wb   =             (xt-xa)*dwb
wb_tl=(-xa_tl)*dwb+(xt-xa)*dwb_tl
IF    (wb > one)THEN
   wb    =one
   dwb   =zero
   wb_tl =zero
   dwb_tl=zero
ELSEIF(wb < zero)THEN
   wb    =zero
   dwb   =zero
   wb_tl =zero
   dwb_tl=zero

ENDIF

bw    =bw    -aw
bw_tl =bw_tl -aw_tl
dbw   =dbw   -daw
dbw_tl=dbw_tl-daw_tl

!w=aw+wb*bw
dw   =daw   + wb   *dbw           + dwb   *bw
dw_tl=daw_tl+(wb_tl*dbw+wb*dbw_tl)+(dwb_tl*bw+dwb*bw_tl)
end subroutine lag_interp_smthweights_tl


!============================================================================
SUBROUTINE lag_interp_smthweights_ad(x,x_AD,xt,aq,aq_AD,bq,bq_AD,w_AD,dw,dw_AD,n)
!============================================================================
INTEGER            ,INTENT(IN   ) :: n
REAL(kind_real)               ,INTENT(IN   ) :: xt
REAL(kind_real)  ,INTENT(IN   ) :: x(n)
REAL(kind_real), INTENT(IN   ) :: aq(n-1),bq(n-1)
REAL(kind_real),INTENT(INOUT) :: aq_ad(n-1),bq_ad(n-1)
REAL(kind_real),INTENT(  OUT) :: dw(n)
REAL(kind_real),INTENT(INOUT) :: x_ad(n),dw_ad(n),w_ad(n)
!-----------------------------------------------------------------------------
REAL(kind_real)                 :: aw(n),bw(n),daw(n),dbw(n)
REAL(kind_real)                 :: aw_ad(n),bw_ad(n),daw_ad(n),dbw_ad(n)
REAL(kind_real)                              :: xa,xb,dwb,wb
REAL(kind_real)                              :: xa_ad,xb_ad,wb_ad,dwb_ad
INTEGER                           :: na
!============================================================================
daw_ad=zero;wb_ad=zero;dbw_ad=zero;dwb_ad=zero;bw_ad=zero
aw_ad =zero;xb_ad=zero;xa_ad =zero

! passive variables needed (from forward code)
CALL lag_interp_weights(x(1:n-1),xt,aq,aw(1:n-1),daw(1:n-1),n-1)
CALL lag_interp_weights(x(2:n  ),xt,bq,bw(2:n  ),dbw(2:n  ),n-1)

aw(n) =zero
daw(n)=zero
bw(1) =zero
dbw(1)=zero
na=n/2
IF(na*2 /= n)stop 'In lag_interp_smthWeights; n must be even'
xa =x(na     )
xb =x(na+1)
dwb=one/(xb-xa)
wb =(xt-xa)*dwb
IF    (wb>one)THEN
   wb =one
   dwb=zero
ELSEIF(wb<zero)THEN
   wb =zero
   dwb=zero
ENDIF
bw =bw -aw
dbw=dbw-daw
!w=aw+wb*bw ! not needed
dw=daw+wb*dbw+dwb*bw ! not needed
!
!
daw_ad=daw_ad+dw_ad
wb_ad =wb_ad +dot_product(dbw,dw_ad)
dbw_ad=dbw_ad+wb*dw_ad
dwb_ad=dwb_ad+dot_product(bw,dw_ad)
bw_ad =bw_ad +dwb*dw_ad

aw_ad=aw_ad+w_ad
wb_ad=wb_ad+dot_product(bw,w_ad)
bw_ad=bw_ad+wb*w_ad

daw_ad=daw_ad-dbw_ad
aw_ad =aw_ad -bw_ad

IF(wb>one)THEN
   wb_ad =zero
   dwb_ad=zero
ELSEIF(wb<zero)THEN
   wb_ad =zero
   dwb_ad=zero
ENDIF

xa_ad        =xa_ad-wb_ad*dwb
dwb_ad       =dwb_ad+(xt-xa)*wb_ad
xb_ad        =xb_ad-(dwb_ad/(xb-xa)**2)
xa_ad        =xa_ad+(dwb_ad/(xb-xa)**2)
x_ad(na+1)=x_ad(na+1)+xb_ad
x_ad(na  )=x_ad(na  )+xa_ad

dbw_ad(1)=zero;bw_ad(1)=zero
daw_ad(n)=zero;aw_ad(n)=zero

CALL lag_interp_weights_ad(x(2:n     ),x_ad(2:n     ),xt,bq,bq_ad,bw_ad(2:n     ),&
             dbw_ad(2:n     ),n-1)

CALL lag_interp_weights_ad(x(1:n-1),x_ad(1:n-1),xt,aq,aq_ad,aw_ad(1:n-1),&
             daw_ad(1:n-1),n-1)

end subroutine lag_interp_smthweights_ad
!============================================================================
end module lag_interp_mod