Computes a single basis function, .
Computation results in triangular tables of the form;
Notice that the position and relative number of nonzero entries in the table depend on "" and on the position of the "" in the first column. This algorithm computes only the nonzero entries.
The value is returned in real :: Nip
variable. The algorithm is similar to Algorithm A2. 2 in its use of the variables temp
and saved
.
Note: This algorithm does not apply to rational basis functions.
Let , , and . The computation of yields
yields
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
class(basis), | intent(in) | :: | me | |||
integer, | intent(in) | :: | i | Basis function index |
||
real(kind=wp), | intent(in) | :: | u | Given value |
pure elemental function OneBasisFun(me,i,u) result(Nip)
class(basis), intent(in) :: me
integer, intent(in) :: i !< Basis function index
real(wp), intent(in) :: u !< Given \( u \) value
real(wp) :: Nip !< \( N_{i,p}(u)\)
!locals
integer :: j, k
real(wp) :: N(0:me%p), Uleft, Uright, saved, temp
!initialize arrays
! if(me%rational)then
! print*, "'OneBasisFun' does not support for rational basis, you already have to compute"
! print*, " all non-zero function for rational basis. Call 'BasisFuns'"
! Nip = 0.0_wp
! return
! end if
N = 0.0_wp
associate( p => me%p, Ui => me%kv, m => me%m ) !get rid of % signs
if (( i == 0 .and. u == Ui(0)) .or. (i == m-p-1 .and. u == Ui(m))) then !/* Special cases */
Nip = 1.0_wp; return
end if
if (u < Ui(i) .or. u >= Ui(i+p+1)) then ! /* Local property */
Nip = 0.0_wp; return
end if
do j=0,p ! /* Initialize zeroth-degree functs */
if (u >= Ui(i+j) .and. u < Ui(i+j+1))then
N(j) = 1.0_wp
else
N(j) = 0.0_wp
end if
end do
do k=1,p ! /* Compute triangular table */
if (N(0) == 0.0_wp)then
saved = 0.0_wp
else
saved = (( u-Ui(i) ) * N(0) ) / ( Ui(i+k) - Ui(i) )
end if
do j=0,p-k
Uleft = Ui(i+j+1)
Uright = Ui(i+j+k+1)
if (N(j+1) == 0.0_wp)then
N(j) = saved
saved = 0.0_wp
else
temp = N(j+1)/(Uright-Uleft)
N(j) = saved+(Uright-u)*temp
saved = (u-Uleft)*temp
end if
end do
end do
Nip = N(0)
end associate
end function OneBasisFun