dlpdp Subroutine

private subroutine dlpdp(a, mda, m, n1, n2, prgopt, x, wnorm, mode, ws, is)

Determine an N1-vector W, and an N2-vector Z which minimizes the Euclidean length of W subject to GW+HZ >= Y. This is the least projected distance problem, LPDP. The matrices G and H are of respective dimensions M by N1 and M by N2.

Called by subprogram DLSI.

  The matrix
             (G H Y)

  occupies rows 1,...,M and cols 1,...,N1+N2+1 of A(*,*).

  The solution (W) is returned in X(*).
               (Z)

Revision history

  • 790701 DATE WRITTEN. Hanson, R. J., (SNLA), Haskell, K. H., (SNLA)
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900328 Added TYPE section. (WRB)
  • 910408 Updated the AUTHOR section. (WRB)

Arguments

Type IntentOptional Attributes Name
real(kind=wp) :: a(mda,*)

A(MDA,N+1), where N=N1+N2.
integer, intent(in) :: mda
integer :: m
integer, intent(in) :: n1
integer, intent(in) :: n2
real(kind=wp) :: prgopt(*)
real(kind=wp) :: x(*)

X(N), where N=N1+N2.
real(kind=wp) :: wnorm
integer, intent(out) :: mode

The value of MODE indicates the status of the computation after returning to the user.

  • MODE=1 The solution was successfully obtained.
  • MODE=2 The inequalities are inconsistent.
real(kind=wp) :: ws(*)

WS((M+2)*(N+7)), where N=N1+N2. This is a slight overestimate for WS(*).
integer :: is(*)

IS(M+N+1), where N=N1+N2.

Calls

proc~~dlpdp~~CallsGraph proc~dlpdp bspline_defc_module::dlpdp proc~dcopy bspline_blas_module::dcopy proc~dlpdp->proc~dcopy proc~ddot bspline_blas_module::ddot proc~dlpdp->proc~ddot proc~dnrm2 bspline_blas_module::dnrm2 proc~dlpdp->proc~dnrm2 proc~dscal bspline_blas_module::dscal proc~dlpdp->proc~dscal proc~dwnnls bspline_defc_module::dwnnls proc~dlpdp->proc~dwnnls proc~dwnlsm bspline_defc_module::dwnlsm proc~dwnnls->proc~dwnlsm proc~dwnlsm->proc~dcopy proc~dwnlsm->proc~dnrm2 proc~dwnlsm->proc~dscal proc~dasum bspline_blas_module::dasum proc~dwnlsm->proc~dasum proc~daxpy bspline_blas_module::daxpy proc~dwnlsm->proc~daxpy proc~dh12 bspline_defc_module::dh12 proc~dwnlsm->proc~dh12 proc~drotm bspline_blas_module::drotm proc~dwnlsm->proc~drotm proc~drotmg bspline_blas_module::drotmg proc~dwnlsm->proc~drotmg proc~dswap bspline_blas_module::dswap proc~dwnlsm->proc~dswap proc~dwnlit bspline_defc_module::dwnlit proc~dwnlsm->proc~dwnlit proc~idamax bspline_blas_module::idamax proc~dwnlsm->proc~idamax proc~dh12->proc~ddot proc~dh12->proc~daxpy proc~dh12->proc~dswap proc~dwnlit->proc~dcopy proc~dwnlit->proc~dscal proc~dwnlit->proc~dh12 proc~dwnlit->proc~drotm proc~dwnlit->proc~drotmg proc~dwnlit->proc~dswap proc~dwnlit->proc~idamax proc~dwnlt1 bspline_defc_module::dwnlt1 proc~dwnlit->proc~dwnlt1 proc~dwnlt2 bspline_defc_module::dwnlt2 proc~dwnlit->proc~dwnlt2 proc~dwnlt3 bspline_defc_module::dwnlt3 proc~dwnlit->proc~dwnlt3 proc~dwnlt1->proc~idamax proc~dwnlt3->proc~dswap

Called by

proc~~dlpdp~~CalledByGraph proc~dlpdp bspline_defc_module::dlpdp proc~dlsi bspline_defc_module::dlsi proc~dlsi->proc~dlpdp proc~dlsei bspline_defc_module::dlsei proc~dlsei->proc~dlsi proc~dfcmn bspline_defc_module::dfcmn proc~dfcmn->proc~dlsei proc~dfc bspline_defc_module::dfc proc~dfc->proc~dfcmn

Source Code

   subroutine dlpdp (a, mda, m, n1, n2, prgopt, x, wnorm, mode, ws, is)

   integer,intent(in) :: mda
   integer :: m
   integer,intent(in) :: n1
   integer,intent(in) :: n2
   real(wp) :: a(mda,*) !! `A(MDA,N+1)`, where `N=N1+N2`.
   real(wp) :: prgopt(*)
   real(wp) :: x(*) !! `X(N)`, where `N=N1+N2`.
   real(wp) :: wnorm
   integer,intent(out) :: mode !! The value of MODE indicates the status of
                               !! the computation after returning to the user.
                               !!
                               !!  * `MODE=1` The solution was successfully obtained.
                               !!  * `MODE=2` The inequalities are inconsistent.
   real(wp) :: ws(*) !! `WS((M+2)*(N+7))`, where `N=N1+N2`. This is a slight overestimate for WS(*).
   integer :: is(*) !! `IS(M+N+1)`, where `N=N1+N2`.

   integer :: i, iw, ix, j, l, modew, n, np1
   real(wp) :: rnorm, sc, ynorm

   real(wp),parameter :: zero = 0.0_wp
   real(wp),parameter :: one = 1.0_wp
   real(wp),parameter :: fac = 0.1_wp

   n = n1 + n2
   mode = 1
   if (m <= 0) then
      if (n > 0) then
         x(1) = zero
         call dcopy(n,x,0,x,1)
      end if
      wnorm = zero
      return
   end if

      np1 = n + 1

      ! SCALE NONZERO ROWS OF INEQUALITY MATRIX TO HAVE LENGTH ONE.
      do i = 1, m
         sc = dnrm2(n,a(i,1),mda)
         if (sc /= zero) then
            sc = one/sc
            call dscal(np1,sc,a(i,1),mda)
         end if
      end do

      ! SCALE RT.-SIDE VECTOR TO HAVE LENGTH ONE (OR ZERO).
      ynorm = dnrm2(m,a(1,np1),1)
      if (ynorm /= zero) then
         sc = one/ynorm
         call dscal(m,sc,a(1,np1),1)
      end if

      ! SCALE COLS OF MATRIX H.
      j = n1 + 1
      do
         if (j > n) exit
         sc = dnrm2(m,a(1,j),1)
         if (sc /= zero) sc = one/sc
         call dscal(m,sc,a(1,j),1)
         x(j) = sc
         j = j + 1
      end do

      if (n1 > 0) then

         ! COPY TRANSPOSE OF (H G Y) TO WORK ARRAY WS(*).
         iw = 0
         do i = 1, m

            ! MOVE COL OF TRANSPOSE OF H INTO WORK ARRAY.
            call dcopy(n2,a(i,n1+1),mda,ws(iw+1),1)
            iw = iw + n2

            ! MOVE COL OF TRANSPOSE OF G INTO WORK ARRAY.
            call dcopy(n1,a(i,1),mda,ws(iw+1),1)
            iw = iw + n1

            ! MOVE COMPONENT OF VECTOR Y INTO WORK ARRAY.
            ws(iw+1) = a(i,np1)
            iw = iw + 1
         end do
         ws(iw+1) = zero
         call dcopy(n,ws(iw+1),0,ws(iw+1),1)
         iw = iw + n
         ws(iw+1) = one
         iw = iw + 1

         ! SOLVE EU=F SUBJECT TO (TRANSPOSE OF H)U=0, U>=0.  THE
         ! MATRIX E = TRANSPOSE OF (G Y), AND THE (N+1)-VECTOR
         ! F = TRANSPOSE OF (0,...,0,1).
         ix = iw + 1
         iw = iw + m

         ! DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF
         ! DWNNLS( ).
         is(1) = 0
         is(2) = 0
         call dwnnls(ws,np1,n2,np1-n2,m,0,prgopt,ws(ix),rnorm, &
                     modew,is,ws(iw+1))

         ! COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY W.
         sc = one - ddot(m,a(1,np1),1,ws(ix),1)
         if (one + fac*abs(sc) == one .or. rnorm <= zero) then
            mode = 2
            return
         end if
         sc = one/sc
         do j = 1, n1
            x(j) = sc*ddot(m,a(1,j),1,ws(ix),1)
         end do
         ! COMPUTE THE VECTOR Q=Y-GW.  OVERWRITE Y WITH THIS
         ! VECTOR.
         do i = 1, m
            a(i,np1) = a(i,np1) - ddot(n1,a(i,1),mda,x,1)
         end do

      end if

      if (n2 > 0) then

         ! COPY TRANSPOSE OF (H Q) TO WORK ARRAY WS(*).
         iw = 0
         do i = 1, m
            call dcopy(n2,a(i,n1+1),mda,ws(iw+1),1)
            iw = iw + n2
            ws(iw+1) = a(i,np1)
            iw = iw + 1
         end do
         ws(iw+1) = zero
         call dcopy(n2,ws(iw+1),0,ws(iw+1),1)
         iw = iw + n2
         ws(iw+1) = one
         iw = iw + 1
         ix = iw + 1
         iw = iw + m

         ! SOLVE RV=S SUBJECT TO V>=0.  THE MATRIX R =(TRANSPOSE
         ! OF (H Q)), WHERE Q=Y-GW.  THE (N2+1)-VECTOR S =(TRANSPOSE
         ! OF (0,...,0,1)).
         !
         ! DO NOT CHECK LENGTHS OF WORK ARRAYS IN THIS USAGE OF
         ! DWNNLS( ).
         is(1) = 0
         is(2) = 0
         call dwnnls(ws,n2+1,0,n2+1,m,0,prgopt,ws(ix),rnorm,modew, &
                     is,ws(iw+1))

         ! COMPUTE THE COMPONENTS OF THE SOLN DENOTED ABOVE BY Z.
         sc = one - ddot(m,a(1,np1),1,ws(ix),1)
         if (one + fac*abs(sc) == one .or. rnorm <= zero) then
            mode = 2
            return
         end if
         sc = one/sc
         do j = 1, n2
            l = n1 + j
            x(l) = sc*ddot(m,a(1,l),1,ws(ix),1)*x(l)
         end do
      end if

      ! ACCOUNT FOR SCALING OF RT.-SIDE VECTOR IN SOLUTION.
      call dscal(n,ynorm,x,1)
      wnorm = dnrm2(n1,x,1)

   end subroutine dlpdp