gjp_algo665.f90

Go to the documentation of this file.

! BEGIN_HEADER
! -----------------------------------------------------------------------------
! Gauss-Jacobi Quadrature Implementation
! Authors: Rohit Goswami <rgoswami[at]ieee.org>
! Source: GaussJacobiQuad Library
! License: MIT
! GitHub Repository: https://github.com/HaoZeke/GaussJacobiQuad
! Date: 2023-08-28
! Commit: c442f77
! -----------------------------------------------------------------------------
! This code is part of the GaussJacobiQuad library, providing an efficient
! implementation for Gauss-Jacobi quadrature nodes and weights computation.
! -----------------------------------------------------------------------------
! To cite this software:
! Rohit Goswami (2023). HaoZeke/GaussJacobiQuad: v0.1.0.
! Zenodo: https://doi.org/10.5281/ZENODO.8285112
! ---------------------------------------------------------------------
! END_HEADER
module gjp_algo665
use gjp_types, only: dp, gjp_sparse_matrix
use gjp_imtqlx, only: imtqlx
use gjp_common, only: jacobi_matrix, jacobi_zeroeth_moment
implicit none
contains
 
subroutine gauss_jacobi_algo665(npts, alpha, beta, x, wts)
    integer, intent(in) :: npts
    real(dp), intent(in) :: alpha, beta
    real(dp), intent(out) :: x(npts), wts(npts)
    real(dp) :: zeroeth_moment
    type(gjp_sparse_matrix) :: jacobi_mat
    real(dp) :: diagonal_elements(npts), &
                off_diagonal_elements(npts - 1)
 
    jacobi_mat = jacobi_matrix(npts, alpha, beta)
    zeroeth_moment = jacobi_zeroeth_moment(alpha, beta)
 
    ! Extract diagonal and off-diagonal elements
    diagonal_elements = jacobi_mat%diagonal(1:npts)
    off_diagonal_elements = jacobi_mat%off_diagonal(1:npts - 1)
 
    ! Initialize weights and knot points
    wts = 0.0_dp
    x = diagonal_elements
    wts(1) = sqrt(zeroeth_moment)
 
    ! Diagonalize the Jacobi matrix using the modified implicitly shifted QL method.
    call imtqlx(npts, x, off_diagonal_elements, wts)
 
    ! The weights are related to the squares of the the eigenvectors
    wts = wts**2
 
end subroutine gauss_jacobi_algo665
 
end module gjp_algo665