lsq Subroutine

private subroutine lsq(m, meq, n, nl, la, l, g, a, b, xl, xu, x, y, w, mode, max_iter_ls, nnls_mode, infbnd)

Minimize with respect to , with upper triangular matrix , and vector , where the unit lower tridiangular matrix is stored columnwise dense in the array with vector stored in its 'diagonal' thus substituting the one-elements of

subject to:

  • ,
  • ,
  • ,

On entry, the user has to provide the arrays l, g, a, b, xl, xu. with dimensions: l(n*(n+1)/2), g(n), a(la,n), b(m), xl(n), xu(n).

The working array w must have at least the following dimension: dim(w) =

  • (3*n+m)*(n+1) for lsq
  • +(n-meq+1)*(mineq+2) + 2*mineq for lsi
  • +(n+mineq)*(n-meq) + 2*meq + n for lsei

with mineq = m - meq + 2*n

On return, no array will be changed by the subroutine.

History

  • coded dieter kraft, april 1987
  • revised march 1989

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: m
integer, intent(in) :: meq
integer, intent(in) :: n
integer, intent(in) :: nl
integer, intent(in) :: la
real(kind=wp), dimension(nl) :: l
real(kind=wp), dimension(n) :: g
real(kind=wp), dimension(la,n) :: a
real(kind=wp), dimension(la) :: b
real(kind=wp), dimension(n) :: xl
real(kind=wp), dimension(n) :: xu
real(kind=wp), dimension(n) :: x

stores the n-dimensional solution vector
real(kind=wp), dimension(m+n+n) :: y

stores the vector of lagrange multipliers of dimension m+n+n (constraints+lower+upper bounds)
real(kind=wp), dimension(*) :: w
integer :: 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
real(kind=wp), intent(in) :: infbnd

"infinity" for the upper and lower bounds.

Calls

proc~~lsq~~CallsGraph proc~lsq slsqp_core::lsq proc~dcopy slsqp_support::dcopy proc~lsq->proc~dcopy proc~ddot slsqp_support::ddot proc~lsq->proc~ddot proc~dscal slsqp_support::dscal proc~lsq->proc~dscal proc~enforce_bounds slsqp_core::enforce_bounds proc~lsq->proc~enforce_bounds proc~lsei slsqp_core::lsei proc~lsq->proc~lsei proc~lsei->proc~dcopy 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~~lsq~~CalledByGraph proc~lsq slsqp_core::lsq 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 lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infBnd)

    implicit none

    integer,intent(in)        :: m
    integer,intent(in)        :: n
    integer,intent(in)        :: meq
    integer,intent(in)        :: nl
    integer,intent(in)        :: la
    real(wp),dimension(n)     :: x  !! stores the n-dimensional solution vector
    real(wp),dimension(m+n+n) :: y  !! stores the vector of lagrange multipliers of dimension
                                    !! m+n+n (constraints+lower+upper bounds)
    integer :: 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
    real(wp),intent(in) :: infbnd !! "infinity" for the upper and lower bounds.

    real(wp),dimension(nl)   :: l
    real(wp),dimension(n)    :: g
    real(wp),dimension(la,n) :: a
    real(wp),dimension(la)   :: b
    real(wp),dimension(*)    :: w
    real(wp),dimension(n)    :: xl
    real(wp),dimension(n)    :: xu
    real(wp) :: diag , xnorm
    integer :: i , ic , id , ie , if , ig , ih , il , im , ip , &
               iu , iw , i1 , i2 , i3 , i4 , mineq , &
               m1 , n1 , n2 , n3, num_unbounded, j

      n1 = n + 1
      mineq = m - meq
      m1 = mineq + n + n

      !  determine whether to solve problem
      !  with inconsistent linerarization (n2=1)
      !  or not (n2=0)

      n2 = n1*n/2 + 1
      if ( n2==nl ) then
         n2 = 0
      else
         n2 = 1
      end if
      n3 = n - n2

      !  recover matrix e and vector f from l and g

      i2 = 1
      i3 = 1
      i4 = 1
      ie = 1
      if = n*n + 1
      do i = 1 , n3
         i1 = n1 - i
         diag = sqrt(l(i2))
         w(i3) = zero
         call dcopy(i1,w(i3),0,w(i3),1)
         call dcopy(i1-n2,l(i2),1,w(i3),n)
         call dscal(i1-n2,diag,w(i3),n)
         w(i3) = diag
         w(if-1+i) = (g(i)-ddot(i-1,w(i4),1,w(if),1))/diag
         i2 = i2 + i1 - n2
         i3 = i3 + n1
         i4 = i4 + n
      end do
      if ( n2==1 ) then
         w(i3) = l(nl)
         w(i4) = zero
         call dcopy(n3,w(i4),0,w(i4),1)
         w(if-1+n) = zero
      end if
      call dscal(n,-one,w(if),1)

      ic = if + n
      id = ic + meq*n

      if ( meq>0 ) then

         !  recover matrix c from upper part of a

         do i = 1 , meq
            call dcopy(n,a(i,1),la,w(ic-1+i),meq)
         end do

         !  recover vector d from upper part of b

         call dcopy(meq,b(1),1,w(id),1)
         call dscal(meq,-one,w(id),1)

      end if

      ig = id + meq

      if ( mineq>0 ) then
         ! recover matrix g from lower part of a
         ! The matrix G(mineq+2*n,m1) is stored at w(ig)
         ! Not all rows will be filled if some of the upper/lower
         ! bounds are unbounded.
         do i = 1 , mineq
            call dcopy(n,a(meq+i,1),la,w(ig-1+i),m1)
         end do
      end if

      ih = ig + m1*n
      iw = ih + mineq + 2*n

      if ( mineq>0 ) then
         ! recover h from lower part of b
         ! The vector H(mineq+2*n) is stored at w(ih)
         call dcopy(mineq,b(meq+1),1,w(ih),1)
         call dscal(mineq,-one,w(ih),1)
      end if

      !  augment matrix g by +i and -i, and,
      !  augment vector h by xl and xu
      !  NaN or infBnd value indicates no bound

      ip = ig + mineq
      il = ih + mineq
      num_unbounded = 0

      do i=1,n
         if (ieee_is_nan(xl(i)) .or. xl(i)<=-infbnd) then
            num_unbounded = num_unbounded + 1
         else
            call update_w(xl(i), one)
         end if
      end do

      do i=1,n
         if (ieee_is_nan(xu(i)) .or. xu(i)>=infbnd) then
            num_unbounded = num_unbounded + 1
         else
            call update_w(xu(i), -one)
         end if
      end do

      call lsei(w(ic),w(id),w(ie),w(if),w(ig),w(ih),max(1,meq),meq,n,n, &
                m1,m1-num_unbounded,n,x,xnorm,w(iw),mode,max_iter_ls,nnls_mode)

      if ( mode==1 ) then
         ! restore lagrange multipliers (only for user-defined variables)
         call dcopy(m,w(iw),1,y(1),1)
         if (n3 > 0) then
            !set rest of the multipliers to nan (they are not used)
            y(m+1) = ieee_value(one, ieee_quiet_nan)
            do i=m+2,m+n3+n3
                y(i) = y(m+1)
            end do
         end if
         call enforce_bounds(x,xl,xu,infbnd)  ! to ensure that bounds are not violated
      end if

      contains

      subroutine update_w(val, fact)
        real(wp),intent(in) :: val  !! xu(i) or xl(i)
        real(wp),intent(in) :: fact !! -1 or 1
        w(il) = fact*val
        do j=1,n
           w(ip + m1*(j-1)) = zero
        end do
        w(ip + m1*(i-1)) = fact
        ip = ip + 1
        il = il + 1
      end subroutine update_w

      end subroutine lsq