function [ x, r ] = dqrlss ( a, lda, m, n, kr, b, jpvt, qraux )

%*****************************************************************************80
%
%% dqrlss() solves a linear system in a least squares sense.
%
%  Discussion:
%
%    DQRLSS must be preceeded by a call to DQRANK.
%
%    The system is to be solved is
%      A * X = B
%    where
%      A is an M by N matrix with rank KR, as determined by DQRANK,
%      B is a given M-vector,
%      X is the N-vector to be computed.
%
%    A solution X, with at most KR nonzero components, is found which
%    minimizes the 2-norm of the residual (A*X-B).
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 August 2016
%
%  Author:
%
%    Original FORTRAN77 version by Dongarra, Moler, Bunch, Stewart.
%    MATLAB version by John Burkardt.
%
%  Input:
%
%    real A(LDA,N), the QR factorization information
%    from DQRANK.  The triangular matrix R of the QR factorization is
%    contained in the upper triangle and information needed to recover
%    the orthogonal matrix Q is stored below the diagonal in A and in
%    the vector QRAUX.
%
%    integer LDA, the leading dimension of A, which must
%    be at least M.
%
%    integer M, the number of rows of A.
%
%    integer N, the number of columns of A.
%
%    integer KR, the rank of the matrix, as estimated
%    by DQRANK.
%
%    real B(M), the right hand side of the linear system.
%
%    integer JPVT(N), the pivot information from DQRANK.
%    Columns JPVT(1), ..., JPVT(KR) of the original matrix are linearly
%    independent to within the tolerance TOL and the remaining columns
%    are linearly dependent.
%
%    real QRAUX(N), auxiliary information from DQRANK
%    defining the QR factorization.
%
%  Output:
%
%    real X(N), a least squares solution to the
%    linear system.
%
%    real R(M), the residual, B - A*X.  RSD may
%    overwite B.
%

%
%  Solve the reduced system of rank KR.
%
  if ( 0 < kr )
    job = 110;
    [ qy, qty, x, r, ab, info ] = dqrsl ( a, lda, m, kr, qraux, b, job );
  end
%
%  Reverse the pivoting to recover the original variable ordering.
%
  jpvt(1:n,1) = - jpvt(1:n,1);

  x(kr+1:n,1) = 0.0;

  for j = 1 : n

    if ( jpvt(j,1) <= 0 )

      k = - jpvt(j,1);
      jpvt(j,1) = k;

      while ( k ~= j )
        t = x(j,1);
        x(j,1) = x(k,1);
        x(k,1) = t;
        jpvt(k,1) = -jpvt(k,1);
        k = jpvt(k,1);
      end

    end

  end

  return
end


