lagzo Subroutine

public subroutine lagzo(n, x, w)

Compute the zeros of Laguerre polynomial Ln(x) in the interval [0,∞], and the corresponding weighting coefficients for Gauss-Laguerre integration

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

Order of the Laguerre polynomial
real(kind=wp), intent(out), dimension(n) :: x

Zeros of the Laguerre polynomial
real(kind=wp), intent(out), dimension(n) :: w

Corresponding weighting coefficients

Source Code

   subroutine lagzo(n,x,w)

      integer,intent(in) :: n !! Order of the Laguerre polynomial
      real(wp),dimension(n),intent(out) :: x !! Zeros of the Laguerre polynomial
      real(wp),dimension(n),intent(out) :: w !! Corresponding weighting coefficients

      real(wp) :: f0 , f1 , fd , gd , hn , p , pd , pf , q , &
                  wp_ , z , z0
      integer :: i , it , j , k , nr

      integer,parameter :: max_iter = 40

      hn = 1.0_wp/n
      pf = 0.0_wp
      pd = 0.0_wp
      do nr = 1 , n
         z = hn
         if ( nr>1 ) z = x(nr-1) + hn*nr**1.27_wp
         it = 0
         do
           it = it + 1
            z0 = z
            p = 1.0_wp
            do i = 1 , nr - 1
               p = p*(z-x(i))
            enddo
            f0 = 1.0_wp
            f1 = 1.0_wp - z
            do k = 2 , n
               pf = ((2.0_wp*k-1.0_wp-z)*f1-(k-1.0_wp)*f0)/k
               pd = k/z*(pf-f1)
               f0 = f1
               f1 = pf
            enddo
            fd = pf/p
            q = 0.0_wp
            do i = 1 , nr - 1
               wp_ = 1.0_wp
               do j = 1 , nr - 1
                  if ( j/=i ) wp_ = wp_*(z-x(j))
               enddo
               q = q + wp_
            enddo
            gd = (pd-q*fd)/p
            z = z - fd/gd
            if ( it<=max_iter .and. abs((z-z0)/z)>1.0e-15_wp ) cycle
            exit
         end do
         x(nr) = z
         w(nr) = 1.0_wp/(z*pd*pd)
      enddo

   end subroutine lagzo