C DQRLSS    SOURCE    PV090527  23/02/07    11:58:37     11592          
      subroutine dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
      
c*****************************************************************************80
c
c! DQRLSS solves a linear system in a least squares sense.
c
c  Discussion:
c
c    DQRLSS must be preceeded by a call to DQRANK.
c
c    The system is to be solved is
c      A * X = B
c    where
c      A is an M by N matrix with rank KR, as determined by DQRANK,
c      B is a given M-vector,
c      X is the N-vector to be computed.
c
c    A solution X, with at most KR nonzero components, is found which
c    minimizes the 2-norm of the residual (A*X-B).
c
c    Once the matrix A has been formed, DQRANK should be
c    called once to decompose it.  Then, for each right hand
c    side B, DQRLSS should be called once to obtain the
c    solution and residual.
c
c  Modified:
c
c    15 April 2012
c
c  Parameters:
c
c    Input, real*8 A(LDA,N), the QR factorization information
c    from DQRANK.  The triangular matrix R of the QR factorization is
c    contained in the upper triangle and information needed to recover
c    the orthogonal matrix Q is stored below the diagonal in A and in
c    the vector QRAUX.
c
c    Input, integer LDA, the leading dimension of A, which must
c    be at least M.
c
c    Input, integer M, the number of rows of A.
c
c    Input, integer N, the number of columns of A.
c
c    Input, integer KR, the rank of the matrix, as estimated
c    by DQRANK.
c
c    Input, real*8 B(M), the right hand side of the linear system.
c
c    Output, real*8 X(N), a least squares solution to the
c    linear system.
c
c    Output, real*8 R(M), the residual, B - A*X.  R may
c    overwite B.
c
c    Input, integer JPVT(N), the pivot information from DQRANK.
c    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
c    independent to within the tolerance TOL and the remaining columns
c    are linearly dependent.
c
c    Input, real*8 QRAUX(N), auxiliary information from DQRANK
c    defining the QR factorization.

        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer lda
        integer m
        integer n
      
        real*8 a(lda,n)
        real*8 b(m)
        integer info
        integer j
        integer job
        integer jpvt(n)
        integer k
        integer kr
        real*8 qraux(n)
        real*8 r(m)
        real*8 t
        real*8 x(n)
      
        if ( kr /= 0 ) then
          job = 110
          call dqrsl ( a, lda, m, kr, qraux, b, r, r, x, r, r, job,info)
        end if
      
        jpvt(1:n) = - jpvt(1:n)
      
        x(kr+1:n) = 0.0D+00
      
        do j = 1, n
      
          if ( jpvt(j) .le. 0 ) then
      
            k = -jpvt(j)
            jpvt(j) = k
      
            do while ( k /= j )
              t = x(j)
              x(j) = x(k)
              x(k) = t
              jpvt(k) = -jpvt(k)
              k = jpvt(k)
            end do
      
          end if
      
        end do
      
        return
      end