lsi Subroutine

private subroutine lsi(e, f, g, h, le, me, lg, mg, n, x, xnorm, w, mode, max_iter_ls, nnls_mode)

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

s.t. .

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

     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) :   n
     dim(w) :   (n+1)*(mg+2) + 2*mg
     dim(jw):   lg

on entry, the user has to provide the arrays 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

  • 03.01.1980, dieter kraft: coded
  • 20.03.1987, dieter kraft: revised to fortran 77

Arguments

Type IntentOptional Attributes Name
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) :: 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) :: xnorm

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

stores the vector of lagrange multipliers in its first 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.
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~~lsi~~CallsGraph proc~lsi slsqp_core::lsi proc~daxpy slsqp_support::daxpy proc~lsi->proc~daxpy proc~ddot slsqp_support::ddot proc~lsi->proc~ddot proc~dnrm2 slsqp_support::dnrm2 proc~lsi->proc~dnrm2 proc~h12 slsqp_core::h12 proc~lsi->proc~h12 proc~ldp slsqp_core::ldp proc~lsi->proc~ldp proc~ldp->proc~daxpy proc~ldp->proc~ddot proc~ldp->proc~dnrm2 proc~bvls_wrapper bvls_module::bvls_wrapper proc~ldp->proc~bvls_wrapper proc~dcopy slsqp_support::dcopy proc~ldp->proc~dcopy 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~~lsi~~CalledByGraph proc~lsi slsqp_core::lsi 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 lsi(e,f,g,h,le,me,lg,mg,n,x,xnorm,w,mode,max_iter_ls,nnls_mode)

    implicit none

    integer,intent(in)                      :: le
    integer,intent(in)                      :: me
    integer,intent(in)                      :: lg
    integer,intent(in)                      :: mg
    integer,intent(in)                      :: n
    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)                    :: xnorm !! stores the residuum of the solution in euclidian norm
    real(wp),dimension(*)   ,intent(inout)  :: w     !! stores the vector of lagrange multipliers in its first
                                                     !! `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.
    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, j
    real(wp) :: t

    !  qr-factors of e and application to f

    do i = 1 , n
        j = min(i+1,n)
        call h12(1,i,i+1,me,e(1,i),1,t,e(1,j),1,le,n-i)
        call h12(2,i,i+1,me,e(1,i),1,t,f,1,1,1)
    end do

    !  transform g and h to get least distance problem

    mode = 5
    do i = 1 , mg
        do j = 1 , n
            if ( abs(e(j,j))<epmach .or. ieee_is_nan(e(j,j))) return
            g(i,j) = (g(i,j)-ddot(j-1,g(i,1),lg,e(1,j),1))/e(j,j)
        end do
        h(i) = h(i) - ddot(n,g(i,1),lg,f,1)
    end do

    !  solve least distance problem

    call ldp(g,lg,mg,n,h,x,xnorm,w,mode,max_iter_ls,nnls_mode)
    if ( mode==1 ) then

        !  solution of original problem

        call daxpy(n,one,f,1,x,1)
        do i = n , 1 , -1
            j = min(i+1,n)
            x(i) = (x(i)-ddot(n-i,e(i,j),le,x(j),1))/e(i,i)
        end do
        j = min(n+1,me)
        t = dnrm2(me-n,f(j),1)
        xnorm = sqrt(xnorm*xnorm+t*t)
    end if

    end subroutine lsi