dlsi Subroutine

   subroutine dlsi (w, mdw, ma, mg, n, prgopt, x, rnorm, mode, ws, ip)

   integer,intent(in) :: mdw !! contain (resp) var. dimension of `W(*,*)`, and matrix dimensions.
   integer,intent(in) :: ma  !! contain (resp) var. dimension of `W(*,*)`, and matrix dimensions.
   integer,intent(in) :: mg  !! contain (resp) var. dimension of `W(*,*)`, and matrix dimensions.
   integer,intent(in) :: n   !! contain (resp) var. dimension of `W(*,*)`, and matrix dimensions.
   real(wp) :: w(mdw,*) !! `W(*,*)` contains:
                        !!
                        !!```
                        !! (A B)
                        !! (G H)
                        !!```
                        !!
                        !! in rows `1,...,MA+MG`,
                        !! cols `1,...,N+1`.
   real(wp),intent(in) :: prgopt(*) !! Program option vector.
   real(wp),intent(out) :: x(*) !! Solution vector(unless MODE=2)
   real(wp),intent(out) :: rnorm !! length of AX-B.
   integer,intent(out) :: mode !! * `=0`   Inequality constraints are compatible.
                               !! * `=2`   Inequality constraints contradictory.
   real(wp) :: ws(*) !! Working storage of dimension `K+N+(MG+2)*(N+7)`,
                     !! where `K=MAX(MA+MG,N)`.
   integer :: ip(*) !! `IP(MG+2*N+1)` Integer working storage

   real(wp) :: anorm, fac, gam, rb, tau, tol, xnorm
   integer :: i, j, k, key, krank, krm1, krp1, l, last, link, m, map1, &
              mdlpdp, minman, n1, n2, n3, next, np1
   logical :: cov, sclcov
   real(wp) :: rnorm_(1) !! JW added for call to [[dhfti]]

   ! Set the nominal tolerance used in the code.
   tol = sqrt(drelpr)

   mode = 0
   rnorm = 0.0_wp
   m = ma + mg
   np1 = n + 1
   krank = 0

   main : block

      if (n<=0 .or. m<=0) exit main

      ! To process option vector.

      cov = .false.
      sclcov = .true.
      last = 1
      link = prgopt(1)

      do
         if (link<=1) exit
         key = prgopt(last+1)
         if (key==1) cov = prgopt(last+2) /= 0.0_wp
         if (key==10) sclcov = prgopt(last+2) == 0.0_wp
         if (key==5) tol = max(drelpr,prgopt(last+2))
         next = prgopt(link)
         last = link
         link = next
      end do

      ! Compute matrix norm of least squares equations.

      anorm = 0.0_wp
      do j = 1,n
         anorm = max(anorm,dasum(ma,w(1,j),1))
      end do

      ! Set tolerance for DHFTI( ) rank test.

      tau = tol*anorm

      ! Compute Householder orthogonal decomposition of matrix.

      call dcopy (n, [0.0_wp], 0, ws, 1)
      call dcopy (ma, w(1, np1), 1, ws, 1)
      k = max(m,n)
      minman = min(ma,n)
      n1 = k + 1
      n2 = n1 + n
      rnorm_(1) = rnorm ! JW
      call dhfti (w, mdw, ma, n, ws, ma, 1, tau, krank, rnorm_, ws(n2), &
                  ws(n1), ip)
      rnorm = rnorm_(1) ! JW
      fac = 1.0_wp
      gam = ma - krank
      if (krank<ma .and. sclcov) fac = rnorm**2/gam

      ! Reduce to DLPDP and solve.

      map1 = ma + 1

      ! Compute inequality rt-hand side for DLPDP.

      if (ma<m) then
         if (minman>0) then
            do i = map1,m
               w(i,np1) = w(i,np1) - ddot(n,w(i,1),mdw,ws,1)
            end do

            ! Apply permutations to col. of inequality constraint matrix.

            do i = 1,minman
               call dswap (mg, w(map1,i), 1, w(map1,ip(i)), 1)
            end do

            ! Apply Householder transformations to constraint matrix.

            if (krank>0 .and. krank<n) then
               do i = krank,1,-1
                  call dh12 (2, i, krank+1, n, w(i,1), mdw, ws(n1+i-1), &
                           w(map1,1), mdw, 1, mg)
               end do
            endif

            ! Compute permuted inequality constraint matrix times r-inv.

            do i = map1,m
               do j = 1,krank
                  w(i,j) = (w(i,j)-ddot(j-1,w(1,j),1,w(i,1),mdw))/w(j,j)
               end do
         end do
         endif

         ! Solve the reduced problem with DLPDP algorithm,
         ! the least projected distance problem.

         call dlpdp(w(map1,1), mdw, mg, krank, n-krank, prgopt, x, &
                  xnorm, mdlpdp, ws(n2), ip(n+1))

         ! Compute solution in original coordinates.

         if (mdlpdp==1) then
            do i = krank,1,-1
               x(i) = (x(i)-ddot(krank-i,w(i,i+1),mdw,x(i+1),1))/w(i,i)
            end do

            ! Apply Householder transformation to solution vector.

            if (krank<n) then
               do i = 1,krank
                  call dh12 (2, i, krank+1, n, w(i,1), mdw, ws(n1+i-1), &
                           x, 1, 1, 1)
               end do
            endif

            ! Repermute variables to their input order.

            if (minman>0) then
               do i = minman,1,-1
                  call dswap (1, x(i), 1, x(ip(i)), 1)
               end do

               ! Variables are now in original coordinates.
               ! Add solution of unconstrained problem.

               do i = 1,n
                  x(i) = x(i) + ws(i)
               end do

               ! Compute the residual vector norm.

               rnorm = sqrt(rnorm**2+xnorm**2)
            endif
         else
            mode = 2
         endif
      else
         call dcopy (n, ws, 1, x, 1)
      endif

      ! Compute covariance matrix based on the orthogonal decomposition
      ! from DHFTI( ).

      if (.not.cov .or. krank<=0) exit main
      krm1 = krank - 1
      krp1 = krank + 1

      ! Copy diagonal terms to working array.

      call dcopy (krank, w, mdw+1, ws(n2), 1)

      ! Reciprocate diagonal terms.

      do j = 1,krank
         w(j,j) = 1.0_wp/w(j,j)
      end do

      ! Invert the upper triangular QR factor on itself.

      if (krank>1) then
         do i = 1,krm1
            do j = i+1,krank
               w(i,j) = -ddot(j-i,w(i,i),mdw,w(i,j),1)*w(j,j)
            end do
         end do
      endif

      ! Compute the inverted factor times its transpose.

      do i = 1,krank
         do j = i,krank
            w(i,j) = ddot(krank+1-j,w(i,j),mdw,w(j,j),mdw)
         end do
      end do

      ! Zero out lower trapezoidal part.
      ! Copy upper triangular to lower triangular part.

      if (krank<n) then
         do j = 1,krank
            call dcopy (j, w(1,j), 1, w(j,1), mdw)
         end do

         do i = krp1,n
            call dcopy (i, [0.0_wp], 0, w(i,1), mdw)
         end do

         ! Apply right side transformations to lower triangle.

         n3 = n2 + krp1
         do i = 1,krank
            l = n1 + i
            k = n2 + i
            rb = ws(l-1)*ws(k-1)

            ! If RB>=0.0_wp, transformation can be regarded as zero.

            if (rb<0.0_wp) then
               rb = 1.0_wp/rb

               ! Store unscaled rank one Householder update in work array.

               call dcopy (n, [0.0_wp], 0, ws(n3), 1)
               l = n1 + i
               k = n3 + i
               ws(k-1) = ws(l-1)

               do j = krp1,n
                  ws(n3+j-1) = w(i,j)
               end do

               do j = 1,n
                  ws(j) = rb*(ddot(j-i,w(j,i),mdw,ws(n3+i-1),1)+ &
                        ddot(n-j+1,w(j,j),1,ws(n3+j-1),1))
               end do

               l = n3 + i
               gam = 0.5_wp*rb*ddot(n-i+1,ws(l-1),1,ws(i),1)
               call daxpy (n-i+1, gam, ws(l-1), 1, ws(i), 1)
               do j = i,n
                  do l = 1,i-1
                     w(j,l) = w(j,l) + ws(n3+j-1)*ws(l)
                  end do
                  do l = i,j
                     w(j,l) = w(j,l) + ws(j)*ws(n3+l-1)+ws(l)*ws(n3+j-1)
                  end do
               end do
            endif
         end do

         ! Copy lower triangle to upper triangle to symmetrize the
         ! covariance matrix.

         do i = 1,n
            call dcopy (i, w(i,1), mdw, w(1,i), 1)
         end do
      endif

      ! Repermute rows and columns.

      do i = minman,1,-1
         k = ip(i)
         if (i/=k) then
            call dswap (1, w(i,i), 1, w(k,k), 1)
            call dswap (i-1, w(1,i), 1, w(1,k), 1)
            call dswap (k-i-1, w(i,i+1), mdw, w(i+1,k), 1)
            call dswap (n-k, w(i, k+1), mdw, w(k, k+1), mdw)
         endif
      end do

      ! Put in normalized residual sum of squares scale factor
      ! and symmetrize the resulting covariance matrix.

      do j = 1,n
         call dscal (j, fac, w(1,j), 1)
         call dcopy (j, w(1,j), 1, w(j,1), mdw)
      end do

   end block main

   ip(1) = krank
   ip(2) = n + max(m,n) + (mg+2)*(n+7)

   end subroutine dlsi