Merge pull request #14 from HaoZeke/varNamesAndDocs

MAINT: Update the variable names and docs

src/GaussJacobiQuad.f90

@@ -24,27 +24,45 @@ module GaussJacobiQuad
 implicit none
 contains
 
-subroutine gauss_jacobi(n, a, b, x, w, method)
-    integer, intent(in) :: n
-    real(dp), intent(in) :: a, b
-    real(dp), intent(out) :: x(n), w(n)
+!> @brief Compute the Gauss-Jacobi quadrature nodes and weights.
+!>
+!> This subroutine calculates the Gauss-Jacobi quadrature nodes and weights for the given parameters @f$\alpha@f$ and @f$\beta@f$,
+!> using the specified method. Gauss-Jacobi quadrature is used to approximate integrals of the form:
+!> \[
+!>   \int_{-1}^{1} (1 - x)^\alpha (1 + x)^\beta f(x) \,dx \approx \sum_{i=1}^{npts} wts_i f(x_i)
+!> \]
+!> where the weights and nodes are calculated with the Jacobi polynomial, which is defined as:
+!> \[
+!>   P_n^{(\alpha, \beta)}(x) = \frac{(\alpha + 1)_n}{n!} \sum_{k=0}^n \binom{n}{k} \frac{(\beta + 1)_{n-k}}{(n-k)!} \left( \frac{x-1}{2} \right)^k \left( \frac{x+1}{2} \right)^{n-k}
+!> \]
+!>
+!> @param[in] npts Number of quadrature points.
+!> @param[in] alpha Parameter alpha in the Jacobi polynomial. Must be greater than -1.
+!> @param[in] beta Parameter beta in the Jacobi polynomial. Must be greater than -1.
+!> @param[out] x Quadrature nodes.
+!> @param[out] wts Quadrature weights.
+!> @param[in] method Method used for calculation. Supported methods are "recurrence" and "gw".
+subroutine gauss_jacobi(npts, alpha, beta, x, wts, method)
+    integer, intent(in) :: npts
+    real(dp), intent(in) :: alpha, beta
+    real(dp), intent(out) :: x(npts), wts(npts)
     character(len=:), allocatable, intent(in) :: method
 
-    if (a <= -1.0_dp) then
+    if (alpha <= -1.0_dp) then
         print*,"Error: alpha must be greater than -1"
         error stop
     end if
 
-    if (b <= -1.0_dp) then
+    if (beta <= -1.0_dp) then
         print*,"Error: beta must be greater than -1"
         error stop
     end if
 
     select case (trim(method))
     case ("recurrence") ! Fails at high beta
-        call gauss_jacobi_rec(n, a, b, x, w)
+        call gauss_jacobi_rec(npts, alpha, beta, x, wts)
     case ("gw") ! Accurate for high beta
-        call gauss_jacobi_gw(n, a, b, x, w)
+        call gauss_jacobi_gw(npts, alpha, beta, x, wts)
     case default
         print*,"Error: Unknown method specified:", method
         print*,"Supported methods: 'recurrence', 'gw''"

src/gjp_gw.f90

@@ -23,7 +23,6 @@
 !>   https://doi.org/10.1007/BF01396453
 module gjp_gw
 use gjp_types, only: dp, gjp_sparse_matrix
-use gjp_constants, only: pi
 use gjp_lapack, only: DSTEQR
 implicit none
 contains
@@ -45,48 +44,52 @@ module gjp_gw
 !> \end{bmatrix}
 !> \]
 !>
-!> `Z` is initialized as the identity matrix, and the eigenvectors are used to compute the weights.
+!> `Z` is initialized as the identity matrix, and the eigenvectors are used to compute the wts.
 !>
-!> @param n Number of nodes
-!> @param a Alpha parameter for Jacobi polynomials
-!> @param b Beta parameter for Jacobi polynomials
-!> @param x (Output) Zeros of Jacobi polynomials
-!> @param w (Output) Weights for Gauss-Jacobi quadrature
-subroutine gauss_jacobi_gw(n, a, b, x, w)
-    integer, intent(in) :: n
-    real(dp), intent(in) :: a, b
-    real(dp), intent(out) :: x(n), w(n)
-    real(dp) :: zmom
-    type(gjp_sparse_matrix) :: jacmat
-    real(dp) :: d(n), e(n - 1), z(n, n), work(2 * n - 2)
-    integer :: info, i
+!> @param[in] npts Number of x
+!> @param[in] alpha parameter for Jacobi polynomials
+!> @param[in] beta parameter for Jacobi polynomials
+!> @param[out] x Zeros of Jacobi polynomials
+!> @param[out] wts weights for Gauss-Jacobi quadrature
+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
 
-    jacmat = jacobi_matrix(n, a, b)
-    zmom = jacobi_zeroeth_moment(a, b)
+    jacobi_mat = jacobi_matrix(npts, alpha, beta)
+    zeroeth_moment = jacobi_zeroeth_moment(alpha, beta)
 
     ! Extract diagonal and off-diagonal elements
-    d = jacmat%diagonal(1:n)
-    e = jacmat%off_diagonal(1:n - 1)
+    diagonal_elements = jacobi_mat%diagonal(1:npts)
+    off_diagonal_elements = jacobi_mat%off_diagonal(1:npts - 1)
 
-    ! Initialize z as identity matrix
-    z = 0.0_dp
-    do i = 1, n
-        z(i, i) = 1.0_dp
+    ! 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.
-    call DSTEQR('V', n, d, e, z, n, work, info)
+    ! Diagonalize the Jacobi matrix using DSTEQR.
+    call DSTEQR('V', npts, diagonal_elements, off_diagonal_elements, &
+                eigenvectors, npts, workspace, computation_info)
 
-    if (info /= 0) then
-        print*,'Error in DSTEQR:', info
-        return
+    if (computation_info /= 0) then
+        write (*, *) 'Error in DSTEQR, info:', computation_info
+        error stop
     end if
 
     ! The eigenvalues are the nodes
-    x = d
+    x = diagonal_elements
     ! The weights are related to the squares of the first components of the
     ! eigenvectors
-    w = z(1, :)**2 * zmom
+    wts = eigenvectors(1, :)**2 * zeroeth_moment
 
 end subroutine gauss_jacobi_gw
 
@@ -101,9 +104,9 @@ end subroutine gauss_jacobi_gw
 !> \end{cases}
 !> \]
 !>
-!> @param n Size of the matrix, number of points
-!> @param alpha Alpha parameter for Jacobi polynomials
-!> @param beta Beta parameter for Jacobi polynomials
+!> @param[in] n Size of the matrix, number of points
+!> @param[in] alpha parameter for Jacobi polynomials
+!> @param[in] beta parameter for Jacobi polynomials
 !> @return A gjp_sparse_matrix representing the Jacobi matrix
 function jacobi_matrix(n, alpha, beta) result(jacmat)
     integer, intent(in) :: n ! Size of the matrix, number of points
@@ -141,8 +144,8 @@ end function jacobi_matrix
 !> \]
 !> Where \(\Gamma\) is the gamma function.
 !>
-!> @param alpha Alpha parameter for Jacobi polynomials
-!> @param beta Beta parameter for Jacobi polynomials
+!> @param[in] alpha parameter for Jacobi polynomials
+!> @param[in] beta parameter for Jacobi polynomials
 !> @return The zeroth moment value
 !> @note The zeroth moment should always be positive
 function jacobi_zeroeth_moment(alpha, beta) result(zmom)

src/gjp_lapack.f90

@@ -35,7 +35,7 @@ module gjp_lapack
 !> @param INFO INFO = 0: successful exit. INFO = K: the algorithm failed to find all the eigenvalues, no eigenvectors or eigenvectors of order K have been found
 interface
     subroutine DSTEQR(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
-        use gjp_types, only: dp
+        import :: dp
         character :: COMPZ
         integer :: N, LDZ, INFO
         real(dp) :: D(*), E(*), Z(LDZ, *), WORK(*)

src/gjp_rec.f90

@@ -28,34 +28,34 @@ module gjp_rec
 contains
 
 ! This returns unsorted roots and weights
-subroutine gauss_jacobi_rec(n, a, b, x, w)
-    integer, intent(in) :: n
-    real(dp), intent(in) :: a, b
-    real(dp), intent(out) :: x(n), w(n)
-    real(dp), dimension(ceiling(n/2._dp)) :: x1, ders1
-    real(dp), dimension(n/2) :: x2, ders2
-    real(dp) :: ders(n), C
+subroutine gauss_jacobi_rec(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), dimension(ceiling(npts/2._dp)) :: x1, ders1
+    real(dp), dimension(npts/2) :: x2, ders2
+    real(dp) :: ders(npts), C
     integer :: idx
-    call recurrence(n, ceiling(n / 2._dp), a, b, x1, ders1)
-    call recurrence(n, n / 2, b, a, x2, ders2)
-    do idx = 1, n / 2
-        x(idx) = -x2(n / 2 - idx + 1)
-        ders(idx) = ders2(n / 2 - idx + 1)
+    call recurrence(npts, ceiling(npts / 2._dp), alpha, beta, x1, ders1)
+    call recurrence(npts, npts / 2, beta, alpha, x2, ders2)
+    do idx = 1, npts / 2
+        x(idx) = -x2(npts / 2 - idx + 1)
+        ders(idx) = ders2(npts / 2 - idx + 1)
     end do
-    do idx = 1, ceiling(n / 2._dp)
-        x(n / 2 + idx) = x1(idx)
-        ders(n / 2 + idx) = ders1(idx)
+    do idx = 1, ceiling(npts / 2._dp)
+        x(npts / 2 + idx) = x1(idx)
+        ders(npts / 2 + idx) = ders1(idx)
     end do
-    w = 1.0d0 / ((1.0d0 - x**2) * ders**2)
-    C = 2**(a + b + 1) * exp(log_gamma(n + a + 1) - &
-                             log_gamma(n + a + b + 1) + &
-                             log_gamma(n + b + 1) - log_gamma(n + 1._dp)); 
-    w = w * C
+    wts = 1.0d0 / ((1.0d0 - x**2) * ders**2)
+    C = 2**(alpha + beta + 1) * exp(log_gamma(npts + alpha + 1) - &
+                                    log_gamma(npts + alpha + beta + 1) + &
+                                    log_gamma(npts + beta + 1) - log_gamma(npts + 1._dp))
+    wts = wts * C
 end subroutine gauss_jacobi_rec
 
-subroutine recurrence(n, n2, a, b, x, PP)
-    integer, intent(in) :: n, n2
-    real(dp), intent(in) :: a, b
+subroutine recurrence(npts, n2, alpha, beta, x, PP)
+    integer, intent(in) :: npts, n2
+    real(dp), intent(in) :: alpha, beta
     real(dp), intent(out) :: x(n2), PP(n2)
     real(dp) :: dx(n2), P(n2)
     integer :: r(n2), l, i
@@ -66,8 +66,8 @@ subroutine recurrence(n, n2, a, b, x, PP)
     end do
 
     do i = 1, n2
-        C = (2 * r(i) + a - 0.5d0) * pi / (2 * n + a + b + 1)
-        T = C + 1 / (2 * n + a + b + 1)**2 * ((0.25d0 - a**2) / tan(0.5d0 * C) - (0.25d0 - b**2) * tan(0.5d0 * C))
+        C = (2 * r(i) + alpha - 0.5d0) * pi / (2 * npts + alpha + beta + 1)
+        T = C + 1 / (2 * npts + alpha + beta + 1)**2 * ((0.25d0 - alpha**2) / tan(0.5d0 * C) - (0.25d0 - beta**2) * tan(0.5d0 * C))
         x(i) = cos(T)
     end do
 
@@ -76,38 +76,38 @@ subroutine recurrence(n, n2, a, b, x, PP)
 
     do while (maxval(abs(dx)) > sqrt(epsilon(1.0d0)) / 1000 .and. l < 10)
         l = l + 1
-        call eval_jacobi_poly(x, n, a, b, P, PP)
+        call eval_jacobi_poly(x, npts, alpha, beta, P, PP)
         dx = -P / PP
         x = x + dx
     end do
 
-    call eval_jacobi_poly(x, n, a, b, P, PP)
+    call eval_jacobi_poly(x, npts, alpha, beta, P, PP)
 end subroutine
 
-subroutine eval_jacobi_poly(x, n, a, b, P, Pp)
-    integer, intent(in) :: n
-    real(dp), intent(in) :: a, b
+subroutine eval_jacobi_poly(x, npts, alpha, beta, P, Pp)
+    integer, intent(in) :: npts
+    real(dp), intent(in) :: alpha, beta
     real(dp), intent(in) :: x(:)
     real(dp), dimension(size(x)), intent(out) :: P, Pp
     real(dp), dimension(size(x)) :: Pm1, Ppm1, Pa1, Ppa1
     integer :: k, i
     real(dp) :: A_val, B_val, C_val, D_val
 
-    P = 0.5d0 * (a - b + (a + b + 2) * x)
+    P = 0.5d0 * (alpha - beta + (alpha + beta + 2) * x)
     Pm1 = 1.0d0
-    Pp = 0.5d0 * (a + b + 2)
+    Pp = 0.5d0 * (alpha + beta + 2)
     Ppm1 = 0.0d0
 
-    if (n == 0) then
+    if (npts == 0) then
         P = Pm1
         Pp = Ppm1
     end if
 
-    do k = 1, n - 1
-        A_val = 2 * (k + 1) * (k + a + b + 1) * (2 * k + a + b)
-        B_val = (2 * k + a + b + 1) * (a**2 - b**2)
-        C_val = product([(2 * k + a + b + i, i=0, 2)])
-        D_val = 2 * (k + a) * (k + b) * (2 * k + a + b + 2)
+    do k = 1, npts - 1
+        A_val = 2 * (k + 1) * (k + alpha + beta + 1) * (2 * k + alpha + beta)
+        B_val = (2 * k + alpha + beta + 1) * (alpha**2 - beta**2)
+        C_val = product([(2 * k + alpha + beta + i, i=0, 2)])
+        D_val = 2 * (k + alpha) * (k + beta) * (2 * k + alpha + beta + 2)
 
         Pa1 = ((B_val + C_val * x) * P - D_val * Pm1) / A_val
         Ppa1 = ((B_val + C_val * x) * Pp + C_val * P - D_val * Ppm1) / A_val