gjp_gw.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_gw
use gjp_types, only: dp, gjp_sparse_matrix
use gjp_lapack, only: dsteqr
use gjp_common, only: jacobi_matrix, jacobi_zeroeth_moment
implicit none
contains
 
subroutine gauss_jacobi_gw(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), &
                eigenvectors(npts, npts), &
                workspace(2 * npts - 2)
    integer :: computation_info, i
 
    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 eigenvectors as identity matrix
    eigenvectors = 0.0_dp
    do i = 1, npts
        eigenvectors(i, i) = 1.0_dp
    end do
 
    ! Diagonalize the Jacobi matrix using DSTEQR.
    call dsteqr('V', npts, diagonal_elements, off_diagonal_elements, &
                eigenvectors, npts, workspace, computation_info)
 
    if (computation_info /= 0) then
        write (*, *) 'Error in DSTEQR, info:', computation_info
        error stop
    end if
 
    ! The eigenvalues are the nodes
    x = diagonal_elements
    ! The weights are related to the squares of the first components of the
    ! eigenvectors
    wts = eigenvectors(1, :)**2 * zeroeth_moment
 
end subroutine gauss_jacobi_gw
 
end module gjp_gw