lsei Subroutine

private subroutine lsei(c, d, e, f, g, h, lc, mc, le, me, lg, mg, n, x, xnrm, w, mode, max_iter_ls, nnls_mode)

for mode=1, the subroutine returns the solution x of equality & inequality constrained least squares problem lsei :

s.t. and .

using QR decomposition & orthogonal basis of nullspace of .

The following dimensions of the arrays defining the problem are necessary:

        dim(c) :   formal (lc,n),    actual (mc,n)
        dim(d) :   formal (lc  ),    actual (mc  )
        dim(e) :   formal (le,n),    actual (me,n)
        dim(f) :   formal (le  ),    actual (me  )
        dim(g) :   formal (lg,n),    actual (mg,n)
        dim(h) :   formal (lg  ),    actual (mg  )
        dim(x) :   formal (n   ),    actual (n   )
        dim(w) :   2*mc+me+(me+mg)*(n-mc)  for lsei
                 +(n-mc+1)*(mg+2)+2*mg     for lsi
        dim(jw):   max(mg,l)

On entry, the user has to provide the arrays C, d, E, f, G, and h. On return, all arrays will be changed by the subroutine.

Reference

  • Chapter 23.6 of Lawson & Hanson: Solving least squares problems.

History

  • 18.5.1981, dieter kraft, dfvlr oberpfaffenhofen
  • 20.3.1987, dieter kraft, dfvlr oberpfaffenhofen

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(inout), dimension(lc,n) :: c
real(kind=wp), intent(inout), dimension(lc) :: d
real(kind=wp), intent(inout), dimension(le,n) :: e
real(kind=wp), intent(inout), dimension(le) :: f
real(kind=wp), intent(inout), dimension(lg,n) :: g
real(kind=wp), intent(inout), dimension(lg) :: h
integer, intent(in) :: lc
integer, intent(in) :: mc
integer, intent(in) :: le
integer, intent(in) :: me
integer, intent(in) :: lg
integer, intent(in) :: mg
integer, intent(in) :: n
real(kind=wp), intent(out), dimension(n) :: x

stores the solution vector
real(kind=wp), intent(out) :: xnrm

stores the residuum of the solution in euclidian norm
real(kind=wp), intent(inout), dimension(*) :: w

on return, stores the vector of lagrange multipliers in its first mc+mg elements
integer, intent(out) :: mode

is a success-failure flag with the following meanings:

  • 1: successful computation,
  • 2: error return because of wrong dimensions (n<1),
  • 3: iteration count exceeded by nnls,
  • 4: inequality constraints incompatible,
  • 5: matrix e is not of full rank,
  • 6: matrix c is not of full rank,
  • 7: rank defect in hfti
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~~lsei~~CallsGraph proc~lsei slsqp_core::lsei proc~dcopy slsqp_support::dcopy proc~lsei->proc~dcopy proc~ddot slsqp_support::ddot proc~lsei->proc~ddot proc~dnrm2 slsqp_support::dnrm2 proc~lsei->proc~dnrm2 proc~h12 slsqp_core::h12 proc~lsei->proc~h12 proc~hfti slsqp_core::hfti proc~lsei->proc~hfti proc~lsi slsqp_core::lsi proc~lsei->proc~lsi proc~hfti->proc~h12 proc~lsi->proc~ddot proc~lsi->proc~dnrm2 proc~lsi->proc~h12 proc~daxpy slsqp_support::daxpy proc~lsi->proc~daxpy proc~ldp slsqp_core::ldp proc~lsi->proc~ldp proc~ldp->proc~dcopy proc~ldp->proc~ddot proc~ldp->proc~dnrm2 proc~ldp->proc~daxpy proc~bvls_wrapper bvls_module::bvls_wrapper proc~ldp->proc~bvls_wrapper proc~nnls slsqp_core::nnls proc~ldp->proc~nnls proc~bvls bvls_module::bvls proc~bvls_wrapper->proc~bvls proc~nnls->proc~h12 proc~g1 slsqp_core::g1 proc~nnls->proc~g1

Called by

proc~~lsei~~CalledByGraph proc~lsei slsqp_core::lsei 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 lsei(c,d,e,f,g,h,lc,mc,le,me,lg,mg,n,x,xnrm,w,mode,&
                    max_iter_ls,nnls_mode)

    implicit none

    integer,intent(in)                      :: lc
    integer,intent(in)                      :: mc
    integer,intent(in)                      :: le
    integer,intent(in)                      :: me
    integer,intent(in)                      :: lg
    integer,intent(in)                      :: mg
    integer,intent(in)                      :: n
    real(wp),dimension(lc,n),intent(inout)  :: c
    real(wp),dimension(lc)  ,intent(inout)  :: d
    real(wp),dimension(le,n),intent(inout)  :: e
    real(wp),dimension(le)  ,intent(inout)  :: f
    real(wp),dimension(lg,n),intent(inout)  :: g
    real(wp),dimension(lg)  ,intent(inout)  :: h
    real(wp),dimension(n)   ,intent(out)    :: x    !! stores the solution vector
    real(wp),intent(out)                    :: xnrm !! stores the residuum of the solution in euclidian norm
    real(wp),dimension(*)   ,intent(inout)  :: w    !! on return, stores the vector of lagrange multipliers
                                                    !! in its first `mc+mg` elements
    integer,intent(out)                     :: mode !! is a success-failure flag with the following meanings:
                                                    !!
                                                    !! * ***1:*** successful computation,
                                                    !! * ***2:*** error return because of wrong dimensions (`n<1`),
                                                    !! * ***3:*** iteration count exceeded by [[nnls]],
                                                    !! * ***4:*** inequality constraints incompatible,
                                                    !! * ***5:*** matrix `e` is not of full rank,
                                                    !! * ***6:*** matrix `c` is not of full rank,
                                                    !! * ***7:*** rank defect in [[hfti]]
    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 , ie, if , ig , iw , j , k , krank , l , mc1
    real(wp) :: t , dum(1)

    mode = 2
    if ( mc<=n ) then
        l = n - mc
        mc1 = mc + 1
        iw = (l+1)*(mg+2) + 2*mg + mc
        ie = iw + mc + 1
        if = ie + me*l
        ig = if + me

        !  triangularize c and apply factors to e and g

        do i = 1 , mc
            j = min(i+1,lc)
            call h12(1,i,i+1,n,c(i,1),lc,w(iw+i),c(j,1),lc,1,mc-i)
            call h12(2,i,i+1,n,c(i,1),lc,w(iw+i),e,le,1,me)
            call h12(2,i,i+1,n,c(i,1),lc,w(iw+i),g,lg,1,mg)
        end do

        !  solve c*x=d and modify f

        mode = 6
        do i = 1 , mc
            if ( abs(c(i,i))<epmach ) return
            x(i) = (d(i)-ddot(i-1,c(i,1),lc,x,1))/c(i,i)
        end do
        mode = 1
        w(mc1) = zero
        !call dcopy(mg-mc,w(mc1),0,w(mc1),1)  ! original code
        call dcopy(mg,w(mc1),0,w(mc1),1)      ! bug fix for when meq = n

        if ( mc/=n ) then

            do i = 1 , me
                w(if-1+i) = f(i) - ddot(mc,e(i,1),le,x,1)
            end do

            !  store transformed e & g

            do i = 1 , me
                call dcopy(l,e(i,mc1),le,w(ie-1+i),me)
            end do
            do i = 1 , mg
                call dcopy(l,g(i,mc1),lg,w(ig-1+i),mg)
            end do

            if ( mg>0 ) then
                !  modify h and solve inequality constrained ls problem

                do i = 1 , mg
                    h(i) = h(i) - ddot(mc,g(i,1),lg,x,1)
                end do
                call lsi(w(ie),w(if),w(ig),h,me,me,mg,mg,l,x(mc1),xnrm,  &
                         w(mc1),mode,max_iter_ls,nnls_mode)

                if ( mc==0 ) return
                t = dnrm2(mc,x,1)
                xnrm = sqrt(xnrm*xnrm+t*t)
                if ( mode/=1 ) return
            else

                ! solve ls without inequality constraints

                mode = 7
                k = max(le,n)
                t = sqrt(epmach)
                call hfti(w(ie),me,me,l,w(if),k,1,t,krank,dum,w,w(l+1))
                xnrm = dum(1)
                call dcopy(l,w(if),1,x(mc1),1)
                if ( krank/=l ) return
                mode = 1
            end if
        end if

        !  solution of original problem and lagrange multipliers

        do i = 1 , me
            f(i) = ddot(n,e(i,1),le,x,1) - f(i)
        end do
        do i = 1 , mc
            d(i) = ddot(me,e(1,i),1,f,1) &
            - ddot(mg,g(1,i),1,w(mc1),1)
        end do

        do i = mc , 1 , -1
            call h12(2,i,i+1,n,c(i,1),lc,w(iw+i),x,1,1,1)
        end do

        do i = mc , 1 , -1
            j = min(i+1,lc)
            w(i) = (d(i)-ddot(mc-i,c(j,i),1,w(j),1))/c(i,i)
        end do
    end if

    end subroutine lsei