slsqp_support.F90 Source File


This file depends on

sourcefile~~slsqp_support.f90~~EfferentGraph sourcefile~slsqp_support.f90 slsqp_support.F90 sourcefile~slsqp_kinds.f90 slsqp_kinds.F90 sourcefile~slsqp_support.f90->sourcefile~slsqp_kinds.f90

Files dependent on this one

sourcefile~~slsqp_support.f90~~AfferentGraph sourcefile~slsqp_support.f90 slsqp_support.F90 sourcefile~bvls_module.f90 bvls_module.f90 sourcefile~bvls_module.f90->sourcefile~slsqp_support.f90 sourcefile~slsqp_core.f90 slsqp_core.f90 sourcefile~slsqp_core.f90->sourcefile~slsqp_support.f90 sourcefile~slsqp_core.f90->sourcefile~bvls_module.f90 sourcefile~slsqp_module.f90 slsqp_module.f90 sourcefile~slsqp_module.f90->sourcefile~slsqp_support.f90 sourcefile~slsqp_module.f90->sourcefile~slsqp_core.f90

Source Code

!*******************************************************************************
!> license: BSD
!
!  Support routines for SLSQP. For example, routines from
!  [BLAS](http://www.netlib.org/blas/) and [LINPACK](http://www.netlib.org/linpack/).
!  These have also been refactored into modern Fortran.

    module slsqp_support

    use slsqp_kinds

    implicit none

    private

    real(wp),parameter,public :: epmach = epsilon(1.0_wp)
    real(wp),parameter,public :: zero   = 0.0_wp
    real(wp),parameter,public :: one    = 1.0_wp
    real(wp),parameter,public :: two    = 2.0_wp
    real(wp),parameter,public :: four   = 4.0_wp
    real(wp),parameter,public :: ten    = 10.0_wp
    real(wp),parameter,public :: hun    = 100.0_wp

    public :: daxpy,dcopy,ddot,dnrm2,dscal

#ifdef HAS_BLAS

    ! linking with an external BLAS library.
    ! Define the interfaces here. Note that this
    ! will only work if the `wp` is the same kind
    ! as used in the BLAS library.

    interface
        pure subroutine daxpy(n,da,dx,incx,dy,incy)
            import :: wp
            implicit none
            integer,intent(in)                  :: n
            real(wp),intent(in)                 :: da
            real(wp),dimension(*),intent(in)    :: dx
            integer,intent(in)                  :: incx
            real(wp),dimension(*),intent(inout) :: dy
            integer,intent(in)                  :: incy
        end subroutine daxpy
        pure subroutine dcopy(n,dx,incx,dy,incy)
            import :: wp
            implicit none
            integer,intent(in)                :: n
            real(wp),dimension(*),intent(in)  :: dx
            integer,intent(in)                :: incx
            real(wp),dimension(*),intent(out) :: dy
            integer,intent(in)                :: incy
        end subroutine dcopy
        pure real(wp) function ddot(n,dx,incx,dy,incy)
            import :: wp
            implicit none
            integer,intent(in)               :: n
            real(wp),dimension(*),intent(in) :: dx
            integer,intent(in)               :: incx
            real(wp),dimension(*),intent(in) :: dy
            integer,intent(in)               :: incy
        end function ddot
        pure function dnrm2(n,x,incx) result(norm)
            import :: wp
            implicit none
            integer,intent(in)               :: incx
            integer,intent(in)               :: n
            real(wp),dimension(*),intent(in) :: x
            real(wp)                         :: norm
        end function dnrm2
        pure subroutine dscal(n,da,dx,incx)
            import :: wp
            implicit none
            integer,intent(in)                  :: n
            real(wp),intent(in)                 :: da
            real(wp),dimension(*),intent(inout) :: dx
            integer,intent(in)                  :: incx
        end subroutine dscal
    end interface

#else

    contains
!*******************************************************************************

!*******************************************************************************
!>
!  constant times a vector plus a vector.
!  uses unrolled loops for increments equal to one.
!
!### Author
!  jack dongarra, linpack, 3/11/78.

    pure subroutine daxpy(n,da,dx,incx,dy,incy)

    implicit none

    integer,intent(in)                  :: n
    real(wp),intent(in)                 :: da
    real(wp),dimension(*),intent(in)    :: dx
    integer,intent(in)                  :: incx
    real(wp),dimension(*),intent(inout) :: dy
    integer,intent(in)                  :: incy

    integer :: i , ix , iy , m , mp1

    if ( n<=0 ) return
    if ( abs(da)<=zero ) return
    if ( incx==1 .and. incy==1 ) then

        ! code for both increments equal to 1

        ! clean-up loop

        m = mod(n,4)
        if ( m/=0 ) then
            do i = 1 , m
                dy(i) = dy(i) + da*dx(i)
            end do
            if ( n<4 ) return
        end if
        mp1 = m + 1
        do i = mp1 , n , 4
            dy(i) = dy(i) + da*dx(i)
            dy(i+1) = dy(i+1) + da*dx(i+1)
            dy(i+2) = dy(i+2) + da*dx(i+2)
            dy(i+3) = dy(i+3) + da*dx(i+3)
        end do

    else

        ! code for unequal increments or equal increments
        ! not equal to 1

        ix = 1
        iy = 1
        if ( incx<0 ) ix = (-n+1)*incx + 1
        if ( incy<0 ) iy = (-n+1)*incy + 1
        do i = 1 , n
            dy(iy) = dy(iy) + da*dx(ix)
            ix = ix + incx
            iy = iy + incy
        end do

    end if

    end subroutine daxpy
!*******************************************************************************

!*******************************************************************************
!>
!  copies a vector, x, to a vector, y.
!  uses unrolled loops for increments equal to one.
!
!### Author
!  jack dongarra, linpack, 3/11/78.

    pure subroutine dcopy(n,dx,incx,dy,incy)

    implicit none

    integer,intent(in)                :: n
    real(wp),dimension(*),intent(in)  :: dx
    integer,intent(in)                :: incx
    real(wp),dimension(*),intent(out) :: dy
    integer,intent(in)                :: incy

    integer :: i , ix , iy , m , mp1

    if ( n<=0 ) return
    if ( incx==1 .and. incy==1 ) then

        ! code for both increments equal to 1

        ! clean-up loop

        m = mod(n,7)
        if ( m/=0 ) then
            do i = 1 , m
                dy(i) = dx(i)
            end do
            if ( n<7 ) return
        end if
        mp1 = m + 1
        do i = mp1 , n , 7
            dy(i) = dx(i)
            dy(i+1) = dx(i+1)
            dy(i+2) = dx(i+2)
            dy(i+3) = dx(i+3)
            dy(i+4) = dx(i+4)
            dy(i+5) = dx(i+5)
            dy(i+6) = dx(i+6)
        end do

    else

        ! code for unequal increments or equal increments
        ! not equal to 1

        ix = 1
        iy = 1
        if ( incx<0 ) ix = (-n+1)*incx + 1
        if ( incy<0 ) iy = (-n+1)*incy + 1
        do i = 1 , n
            dy(iy) = dx(ix)
            ix = ix + incx
            iy = iy + incy
        end do

    end if

    end subroutine dcopy
!*******************************************************************************

!*******************************************************************************
!>
!  forms the dot product of two vectors.
!  uses unrolled loops for increments equal to one.
!
!### Author
!  jack dongarra, linpack, 3/11/78.

    pure real(wp) function ddot(n,dx,incx,dy,incy)

    implicit none

    integer,intent(in)               :: n
    real(wp),dimension(*),intent(in) :: dx
    integer,intent(in)               :: incx
    real(wp),dimension(*),intent(in) :: dy
    integer,intent(in)               :: incy

    real(wp) :: dtemp
    integer :: i , ix , iy , m , mp1

    ddot = zero
    dtemp = zero
    if ( n<=0 ) return
    if ( incx==1 .and. incy==1 ) then

        ! code for both increments equal to 1

        ! clean-up loop

        m = mod(n,5)
        if ( m/=0 ) then
            do i = 1 , m
                dtemp = dtemp + dx(i)*dy(i)
            end do
            if ( n<5 ) then
                ddot = dtemp
                return
            end if
        end if
        mp1 = m + 1
        do i = mp1 , n , 5
            dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + &
                    dx(i+2)*dy(i+2) + dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4)
        end do
        ddot = dtemp

    else

        ! code for unequal increments or equal increments
        ! not equal to 1

        ix = 1
        iy = 1
        if ( incx<0 ) ix = (-n+1)*incx + 1
        if ( incy<0 ) iy = (-n+1)*incy + 1
        do i = 1 , n
            dtemp = dtemp + dx(ix)*dy(iy)
            ix = ix + incx
            iy = iy + incy
        end do
        ddot = dtemp

    end if

    end function ddot
!*******************************************************************************

!*******************************************************************************
!>
!  Function that returns the Euclidean norm
!  \( \sqrt{ \mathbf{x}^T \mathbf{x} } \) of a vector \( \mathbf{x} \).
!
!### Further details
!
!  * this version written on 25-october-1982.
!  * modified on 14-october-1993 to inline the call to dlassq.
!    sven hammarling, nag ltd.
!  * Converted to modern Fortran, Jacob Williams, Jan. 2016.
!
!@note Replaced original SLSQP routine with this one from
!      [BLAS](http://netlib.sandia.gov/blas/dnrm2.f).

    pure function dnrm2(n,x,incx) result(norm)

    implicit none

    integer,intent(in)               :: incx
    integer,intent(in)               :: n
    real(wp),dimension(*),intent(in) :: x
    real(wp)                         :: norm

    real(wp) :: absxi , scale , ssq
    integer :: ix

    if ( n<1 .or. incx<1 ) then
        norm = zero
    elseif ( n==1 ) then
        norm = abs(x(1))
    else
        scale = zero
        ssq = one
        ! the following loop is equivalent to this call to the lapack
        ! auxiliary routine:
        ! call dlassq( n, x, incx, scale, ssq )
        do ix = 1 , 1 + (n-1)*incx , incx
            if ( abs(x(ix))>zero ) then
                absxi = abs(x(ix))
                if ( scale<absxi ) then
                    ssq = one + ssq*(scale/absxi)**2
                    scale = absxi
                else
                    ssq = ssq + (absxi/scale)**2
                end if
            end if
        end do
        norm = scale*sqrt(ssq)
    end if

    end function dnrm2
!*******************************************************************************

!*******************************************************************************
!>
!  scales a vector by a constant.
!  uses unrolled loops for increment equal to one.
!
!### Author
!  jack dongarra, linpack, 3/11/78.

    pure subroutine dscal(n,da,dx,incx)

    implicit none

    integer,intent(in)                  :: n
    real(wp),intent(in)                 :: da
    real(wp),dimension(*),intent(inout) :: dx
    integer,intent(in)                  :: incx

    integer :: i , m , mp1 , nincx

    if ( n<=0 .or. incx<=0 ) return
    if ( incx==1 ) then

        ! code for increment equal to 1

        ! clean-up loop

        m = mod(n,5)
        if ( m/=0 ) then
            do i = 1 , m
                dx(i) = da*dx(i)
            end do
            if ( n<5 ) return
        end if
        mp1 = m + 1
        do i = mp1 , n , 5
            dx(i) = da*dx(i)
            dx(i+1) = da*dx(i+1)
            dx(i+2) = da*dx(i+2)
            dx(i+3) = da*dx(i+3)
            dx(i+4) = da*dx(i+4)
        end do

    else

        ! code for increment not equal to 1

        nincx = n*incx
        do i = 1 , nincx , incx
            dx(i) = da*dx(i)
        end do

    end if

    end subroutine dscal
!*******************************************************************************

#endif

!*******************************************************************************
    end module slsqp_support
!*******************************************************************************