ldp Subroutine

private subroutine ldp(g, mg, m, n, h, x, xnorm, w, mode, max_iter_ls, nnls_mode)

Least distance programming routine. Minimize subject to .

The declared dimension of w must be at least (n+1)*(m+2)+2*m:

       first (n+1)*m locs of w = matrix e for problem nnls.
       next      n+1 locs of w = vector f for problem nnls.
       next      n+1 locs of w = vector z for problem nnls.
       next        m locs of w = vector y for problem nnls.
       next        m locs of w = vector wdual for problem nnls.

References

  • C.L. Lawson, R.J. Hanson, 'Solving least squares problems' Prentice Hall, 1974. (revised 1995 edition)
  • lawson-hanson from Netlib.

History

  • Jacob Williams, refactored into modern Fortran, Jan. 2016.

Note

The 1995 version of this routine may have some sort of problem. Using a refactored version of the original routine.

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in), dimension(mg,n) :: g

on entry g stores the m by n matrix of linear inequality constraints. g has first dimensioning parameter mg
integer, intent(in) :: mg
integer, intent(in) :: m
integer, intent(in) :: n
real(kind=wp), intent(in), dimension(m) :: h

the right side of the inequality system.
real(kind=wp), intent(out), dimension(n) :: x

solution vector x if mode=1.
real(kind=wp), intent(out) :: xnorm

euclidian norm of the solution vector if computation is successful
real(kind=wp), intent(inout), dimension(*) :: w

w is a one dimensional working space, the length of which should be at least (m+2)*(n+1) + 2*m. on exit w stores the lagrange multipliers associated with the constraints. at the solution of problem ldp.
integer, intent(out) :: mode

success-failure flag with the following meanings:

  • 1: successful computation,
  • 2: error return because of wrong dimensions (n<=0),
  • 3: iteration count exceeded by nnls,
  • 4: inequality constraints incompatible.
integer, intent(in) :: max_iter_ls

maximum number of iterations in the nnls problem
integer, intent(in) :: nnls_mode

which NNLS method to use

Calls

proc~~ldp~~CallsGraph proc~ldp slsqp_core::ldp proc~bvls_wrapper bvls_module::bvls_wrapper proc~ldp->proc~bvls_wrapper proc~daxpy slsqp_support::daxpy proc~ldp->proc~daxpy proc~dcopy slsqp_support::dcopy proc~ldp->proc~dcopy proc~ddot slsqp_support::ddot proc~ldp->proc~ddot proc~dnrm2 slsqp_support::dnrm2 proc~ldp->proc~dnrm2 proc~nnls slsqp_core::nnls proc~ldp->proc~nnls proc~bvls bvls_module::bvls proc~bvls_wrapper->proc~bvls proc~g1 slsqp_core::g1 proc~nnls->proc~g1 proc~h12 slsqp_core::h12 proc~nnls->proc~h12

Called by

proc~~ldp~~CalledByGraph proc~ldp slsqp_core::ldp proc~lsi slsqp_core::lsi proc~lsi->proc~ldp proc~lsei slsqp_core::lsei proc~lsei->proc~lsi proc~lsq slsqp_core::lsq proc~lsq->proc~lsei proc~slsqpb slsqp_core::slsqpb proc~slsqpb->proc~lsq proc~slsqp slsqp_core::slsqp proc~slsqp->proc~slsqpb proc~slsqp_wrapper slsqp_module::slsqp_solver%slsqp_wrapper proc~slsqp_wrapper->proc~slsqp

Source Code

    subroutine ldp(g,mg,m,n,h,x,xnorm,w,mode,max_iter_ls,nnls_mode)

    implicit none

    integer,intent(in)                   :: mg
    integer,intent(in)                   :: m
    integer,intent(in)                   :: n
    real(wp),dimension(mg,n),intent(in)  :: g     !! on entry `g` stores the `m` by `n` matrix of
                                                  !! linear inequality constraints. `g` has first
                                                  !! dimensioning parameter `mg`
    real(wp),dimension(m),intent(in)     :: h     !! the right side of the inequality system.
    real(wp),dimension(n),intent(out)    :: x     !! solution vector `x` if `mode=1`.
    real(wp),dimension(*),intent(inout)  :: w     !! `w` is a one dimensional working space, the length
                                                  !! of which should be at least `(m+2)*(n+1) + 2*m`.
                                                  !! on exit `w` stores the lagrange multipliers
                                                  !! associated with the constraints.
                                                  !! at the solution of problem `ldp`.
    real(wp),intent(out)                 :: xnorm !! euclidian norm of the solution vector
                                                  !! if computation is successful
    integer,intent(out)                  :: mode  !! success-failure flag with the following meanings:
                                                  !!
                                                  !! * ***1:*** successful computation,
                                                  !! * ***2:*** error return because of wrong dimensions (`n<=0`),
                                                  !! * ***3:*** iteration count exceeded by [[nnls]],
                                                  !! * ***4:*** inequality constraints incompatible.
    integer,intent(in) :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem
    integer,intent(in) :: nnls_mode !! which NNLS method to use

    integer :: i , iw , iwdual , iy , iz , j , jf , n1
    real(wp) :: fac , rnorm

    if ( n<=0 ) then
       ! error return.
       mode = 2
    else
        ! state dual problem
        mode = 1
        x = zero
        xnorm = zero
        if (m/=0) then
            iw=0
            do j=1,m
                do i=1,n
                    iw=iw+1
                    w(iw)=g(j,i)
                end do
                iw=iw+1
                w(iw)=h(j)
            end do
            jf=iw+1
            do i=1,n
                iw=iw+1
                w(iw)=zero
            end do
            w(iw+1)=one
            n1=n+1
            iz=iw+2
            iy=iz+n1
            iwdual=iy+m
            ! solve dual problem
            select case (nnls_mode)
            case(1) ! original
                call nnls (w,n1,n1,m,w(jf),w(iy),rnorm,w(iwdual),w(iz),mode,max_iter_ls)
            case(2) ! new version
                call bvls_wrapper(w,n1,n1,m,w(jf),w(iy),rnorm,w(iwdual),w(iz),mode,max_iter_ls)
            case default
                error stop 'invalid nnls_mode'
            end select

            if (mode==1) then
                mode=4
                if (rnorm>zero) then
                    !  compute solution of primal problem
                    fac=one-ddot(m,h,1,w(iy),1)
                    if (ieee_is_nan(fac)) return
                    if (fac>=eps) then
                        mode=1
                        fac=one/fac
                        do j=1,n
                            x(j)=fac*ddot(m,g(1,j),1,w(iy),1)
                        end do
                        xnorm=dnrm2(n,x,1)
                        ! compute lagrange multipliers for primal problem
                        w(1)=zero
                        call dcopy(m,w(1),0,w,1)
                        call daxpy(m,fac,w(iy),1,w,1)
                    end if
                end if
            end if
        end if
    end if

    end subroutine ldp