OneBasisFun Function

private pure elemental function OneBasisFun(me, i, u) result(Nip)

Computes a single basis function, .

Implementation Details of Algorithm A2.4

Computation results in triangular tables of the form;

Algorithm 2.4

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.


Ex2.5 Computing One Basis Function

Let , , and . The computation of yields

yields


Arguments

Type IntentOptional AttributesName
class(basis), intent(in) :: me
integer, intent(in) :: i

Basis function index

real(kind=wp), intent(in) :: u

Given value

Return Value real(kind=wp)


Calls

proc~~onebasisfun~~CallsGraph proc~onebasisfun OneBasisFun ui ui proc~onebasisfun->ui

Contents

Source Code


Source Code

        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