C QRSOLV    SOURCE    PV090527  23/01/30    21:15:04     11574          
      subroutine daxpy ( n, da, dx, incx, dy, incy )
      
c*****************************************************************************80
c
c! DAXPY computes constant times a vector plus a vector.
c
c  Discussion:
c
c    This routine uses double precision real arithmetic.
c
c    This routine uses unrolled loops for increments equal to one.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    16 May 2005
c
c  Author:
c
c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
c    David Kincaid, Fred Krogh.
c    FORTRAN90 version by John Burkardt.
c
c  Reference:
c
c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
c    Algorithm 539, 
c    Basic Linear Algebra Subprograms for Fortran Usage,
c    ACM Transactions on Mathematical Software, 
c    Volume 5, Number 3, September 1979, pages 308-323.
c
c  Parameters:
c
c    Input, integer N, the number of elements in DX and DY.
c
c    Input, real*8 DA, the multiplier of DX.
c
c    Input, real*8 DX(*), the first vector.
c
c    Input, integer INCX, the increment between successive 
c    entries of DX.
c
c    Input/output, real*8 DY(*), the second vector.
c    On output, DY(*) has been replaced by DY(*) + DA * DX(*).
c
c    Input, integer INCY, the increment between successive 
c    entries of DY.

        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        real*8 da
        real*8 dx(*)
        real*8 dy(*)
        integer i
        integer incx
        integer incy
        integer ix
        integer iy
        integer m
        integer n
      
        if ( n .le. 0 ) then
          return
        end if
      
        if ( da .eq. 0.0D+00 ) then
          return
        end if

c  Code for unequal increments or equal increments
c  not equal to 1.

        if ( incx /= 1 .or. incy /= 1 ) then
      
          if ( 0 .le. incx ) then
            ix = 1
          else
            ix = ( - n + 1 ) * incx + 1
          end if
      
          if ( 0 .le. incy ) then
            iy = 1
          else
            iy = ( - n + 1 ) * incy + 1
          end if
      
          do i = 1, n
            dy(iy) = dy(iy) + da * dx(ix)
            ix = ix + incx
            iy = iy + incy
          end do

c  Code for both increments equal to 1.

        else
      
          m = mod ( n, 4 )
      
          dy(1:m) = dy(1:m) + da * dx(1:m)
      
          do i = m+1, n, 4
            dy(i  ) = dy(i  ) + da * dx(i  )
            dy(i+1) = dy(i+1) + da * dx(i+1)
            dy(i+2) = dy(i+2) + da * dx(i+2)
            dy(i+3) = dy(i+3) + da * dx(i+3)
          end do
      
        end if
      
        return
      end
      function dnrm2 ( n, x, incx )
      
c*****************************************************************************80
c
c! DNRM2 returns the euclidean norm of a vector.
c
c  Discussion:
c
c    This routine uses double precision real arithmetic.
c
c     DNRM2 ( X ) = sqrt ( X' * X )
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    16 May 2005
c
c  Author:
c
c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
c    David Kincaid, Fred Krogh.
c    FORTRAN90 version by John Burkardt.
c
c  Reference:
c
c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
c    Algorithm 539, 
c    Basic Linear Algebra Subprograms for Fortran Usage,
c    ACM Transactions on Mathematical Software,
c    Volume 5, Number 3, September 1979, pages 308-323.
c
c  Parameters:
c
c    Input, integer N, the number of entries in the vector.
c
c    Input, real*8 X(*), the vector whose norm is to be computed.
c
c    Input, integer INCX, the increment between successive 
c    entries of X.
c
c    Output, real*8 DNRM2, the Euclidean norm of X.
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        real*8 absxi
        real*8 dnrm2
        integer incx
        integer ix
        integer n
        real*8 norm
        real*8 scale
        real*8 ssq
        real*8 x(*)
      
        if ( n .lt. 1 .or. incx .lt. 1 ) then
      
          norm  = 0.0D+00
      
        else if ( n .eq. 1 ) then
      
          norm  = abs ( x(1) )
      
        else
      
          scale = 0.0D+00
          ssq = 1.0D+00
      
          do ix = 1, 1 + ( n - 1 ) * incx, incx
            if ( x(ix) /= 0.0D+00 ) then
              absxi = abs ( x(ix) )
              if ( scale .lt. absxi ) then
                ssq = 1.0D+00 + ssq * ( scale / absxi )**2
                scale = absxi
              else
                ssq = ssq + ( absxi / scale )**2
              end if
            end if
          end do
          norm  = scale * sqrt ( ssq )
        end if
      
        dnrm2 = norm
      
        return
      end
      subroutine dqrank ( a, lda, m, n, tol, kr, jpvt, qraux, work )
      
c*****************************************************************************80
c
c! DQRANK computes the QR factorization of a rectangular matrix.
c
c  Discussion:
c
c    This routine is used in conjunction with sqrlss to solve
c    overdetermined, underdetermined and singular linear systems
c    in a least squares sense.
c
c    DQRANK uses the LINPACK subroutine DQRDC to compute the QR
c    factorization, with column pivoting, of an M by N matrix A.
c    The numerical rank is determined using the tolerance TOL.
c
c    Note that on output, ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
c    of the condition number of the matrix of independent columns,
c    and of R.  This estimate will be .le. 1/TOL.
c
c  Reference:
c
c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c  Parameters:
c
c    Input/output, real*8 A(LDA,N).  On input, the matrix whose
c    decomposition is to be computed.  On output, the information from DQRDC.
c    The triangular matrix R of the QR factorization is contained in the
c    upper triangle and information needed to recover the orthogonal
c    matrix Q is stored below the diagonal in A and in 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, real*8 TOL, a relative tolerance used to determine the
c    numerical rank.  The problem should be scaled so that all the elements
c    of A have roughly the same absolute accuracy, EPS.  Then a reasonable
c    value for TOL is roughly EPS divided by the magnitude of the largest
c    element.
c
c    Output, integer KR, the numerical rank.
c
c    Output, integer JPVT(N), the pivot information from DQRDC.
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    Output, real*8 QRAUX(N), will contain extra information defining
c    the QR factorization.
c
c    Workspace, real*8 WORK(N).
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer lda
        integer n
      
        real*8 a(lda,n)
        integer j
        integer job
        integer jpvt(n)
        integer k
        integer kr
        integer m
        real*8 qraux(n)
        real*8 tol
        real*8 work(n)
      
        jpvt(1:n) = 0
        job = 1
      
        call dqrdc ( a, lda, m, n, qraux, jpvt, work, job )
      
        kr = 0
        k = min ( m, n )
      
        do j = 1, k
          if ( abs ( a(j,j) ) .le. tol * abs ( a(1,1) ) ) then
            return
          end if
          kr = j
        end do
      
        return
      end
      subroutine dqrdc ( a, lda, n, p, qraux, jpvt, work, job )
      
c*****************************************************************************80
c
c! DQRDC computes the QR factorization of a real rectangular matrix.
c
c  Discussion:
c
c    DQRDC uses Householder transformations.
c
c    Column pivoting based on the 2-norms of the reduced columns may be
c    performed at the user's option.
c
c  Reference:
c
c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c  Parameters:
c
c    Input/output, real*8 A(LDA,P).  On input, the N by P matrix
c    whose decomposition is to be computed.  On output, A contains in
c    its upper triangle the upper triangular matrix R of the QR
c    factorization.  Below its diagonal A contains information from
c    which the orthogonal part of the decomposition can be recovered.
c    Note that if pivoting has been requested, the decomposition is not that
c    of the original matrix A but that of A with its columns permuted
c    as described by JPVT.
c
c    Input, integer LDA, the leading dimension of the array A.
c    LDA must be at least N.
c
c    Input, integer N, the number of rows of the matrix A.
c
c    Input, integer P, the number of columns of the matrix A.
c
c    Output, real*8 QRAUX(P), contains further information required
c    to recover the orthogonal part of the decomposition.
c
c    Input/output, integer JPVT(P).  On input, JPVT contains
c    integers that control the selection of the pivot columns.  The K-th
c    column A(*,K) of A is placed in one of three classes according to the
c    value of JPVT(K).
c      .gt. 0, then A(K) is an initial column.
c      = 0, then A(K) is a free column.
c      .lt. 0, then A(K) is a final column.
c    Before the decomposition is computed, initial columns are moved to
c    the beginning of the array A and final columns to the end.  Both
c    initial and final columns are frozen in place during the computation
c    and only free columns are moved.  At the K-th stage of the
c    reduction, if A(*,K) is occupied by a free column it is interchanged
c    with the free column of largest reduced norm.  JPVT is not referenced
c    if JOB .eq. 0.  On output, JPVT(K) contains the index of the column of the
c    original matrix that has been interchanged into the K-th column, if
c    pivoting was requested.
c
c    Workspace, real*8 WORK(P).  WORK is not referenced if JOB .eq. 0.
c
c    Input, integer JOB, initiates column pivoting.
c    0, no pivoting is done.
c    nonzero, pivoting is done.
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer lda
        integer n
        integer p
      
        real*8 a(lda,p)
        integer jpvt(p)
        real*8 qraux(p)
        real*8 work(p)
        integer j
        integer job
        integer jp
        integer l
        integer lup
        integer maxj
        real*8 maxnrm
        real*8 nrmxl
        integer pl
        integer pu
        real*8 ddot
        real*8 dnrm2
        logical swapj
        real*8 t
        real*8 tt
      
        pl = 1
        pu = 0

c  If pivoting is requested, rearrange the columns.

        if ( job /= 0 ) then
      
          do j = 1, p
      
            swapj = 0 .lt. jpvt(j)
      
            if ( jpvt(j) .lt. 0 ) then
              jpvt(j) = - j
            else
              jpvt(j) = j
            end if
      
            if ( swapj ) then
      
              if ( j /= pl ) then
                call dswap ( n, a(1,pl), 1, a(1,j), 1 )
              end if
      
              jpvt(j) = jpvt(pl)
              jpvt(pl) = j
              pl = pl + 1
      
            end if
      
          end do
      
          pu = p
      
          do j = p, 1, -1
      
            if ( jpvt(j) .lt. 0 ) then
      
              jpvt(j) = - jpvt(j)
      
              if ( j /= pu ) then
                call dswap ( n, a(1,pu), 1, a(1,j), 1 )
                jp = jpvt(pu)
                jpvt(pu) = jpvt(j)
                jpvt(j) = jp
              end if
      
              pu = pu - 1
      
            end if
      
          end do
      
        end if

c  Compute the norms of the free columns.

        do j = pl, pu
          qraux(j) = dnrm2 ( n, a(1,j), 1 )
        end do
      
        work(pl:pu) = qraux(pl:pu)

c  Perform the Householder reduction of A.

        lup = min ( n, p )
      
        do l = 1, lup

c  Bring the column of largest norm into the pivot position.

          if ( (pl .le. l) .and. (l .lt. pu) ) then
      
            maxnrm = 0.0D+00
            maxj = l
            do j = l, pu
              if ( maxnrm .lt. qraux(j) ) then
                maxnrm = qraux(j)
                maxj = j
              end if
            end do
      
            if ( maxj /= l ) then
              call dswap ( n, a(1,l), 1, a(1,maxj), 1 )
              qraux(maxj) = qraux(l)
              work(maxj) = work(l)
              jp = jpvt(maxj)
              jpvt(maxj) = jpvt(l)
              jpvt(l) = jp
            end if
      
          end if

c  Compute the Householder transformation for column L.

          qraux(l) = 0.0D+00
      
          if ( l /= n ) then
      
            nrmxl = dnrm2 ( n-l+1, a(l,l), 1 )
      
            if ( nrmxl /= 0.0D+00 ) then
      
              if ( a(l,l) /= 0.0D+00 ) then
                nrmxl = sign ( nrmxl, a(l,l) )
              end if
      
              call dscal ( n-l+1, 1.0D+00 / nrmxl, a(l,l), 1 )
              a(l,l) = 1.0D+00 + a(l,l)

c  Apply the transformation to the remaining columns, updating the norms.

              do j = l + 1, p
      
                t = - ddot ( n-l+1, a(l,l), 1, a(l,j), 1 ) / a(l,l)
                call daxpy ( n-l+1, t, a(l,l), 1, a(l,j), 1 )
      
                if ( (pl .le. j) .and. (j .le. pu) ) then
      
                  if ( qraux(j) /= 0.0D+00 ) then
      
                    tt = 1.0D+00 - ( abs ( a(l,j) ) / qraux(j) )**2
                    tt = max ( tt, 0.0D+00 )
                    t = tt
                    tt = 1.0D+00 + 0.05D+00 * tt * (qraux(j)/work(j))**2
      
                    if ( tt /= 1.0D+00 ) then
                      qraux(j) = qraux(j) * sqrt ( t )
                    else
                      qraux(j) = dnrm2 ( n-l, a(l+1,j), 1 )
                      work(j) = qraux(j)
                    end if
      
                  end if
      
                end if
      
              end do

c  Save the transformation.

              qraux(l) = a(l,l)
              a(l,l) = - nrmxl
      
            end if
      
          end if
      
        end do
      
        return
      end
      subroutine dqrls ( a, lda, m, n, tol, kr, b, x,r,jpvt,qraux,work,
     #  itask, ind )
      
c*****************************************************************************80
c
c! DQRLS factors and solves a linear system in the least squares sense.
c
c  Discussion:
c
c    The linear system may be overdetermined, underdetermined or singular.
c    The solution is obtained using a QR factorization of the
c    coefficient matrix.
c
c    DQRLS can be efficiently used to solve several least squares
c    problems with the same matrix A.  The first system is solved
c    with ITASK = 1.  The subsequent systems are solved with
c    ITASK = 2, to avoid the recomputation of the matrix factors.
c    The parameters KR, JPVT, and QRAUX must not be modified
c    between calls to DQRLS.
c
c    DQRLS is used to solve in a least squares sense
c    overdetermined, underdetermined and singular linear systems.
c    The system is A*X approximates B where A is M by N.
c    B is a given M-vector, and X is the N-vector to be computed.
c    A solution X is found which minimimzes the sum of squares (2-norm)
c    of the residual,  A*X - B.
c
c    The numerical rank of A is determined using the tolerance TOL.
c
c    DQRLS uses the LINPACK subroutine DQRDC to compute the QR
c    factorization, with column pivoting, of an M by N matrix A.
c
c  Modified:
c
c    15 April 2012
c
c  Reference:
c
c    David Kahaner, Cleve Moler, Steven Nash,
c    Numerical Methods and Software,
c    Prentice Hall, 1989,
c    ISBN: 0-13-627258-4,
c    LC: TA345.K34.
c
c  Parameters:
c
c    Input/output, real*8 A(LDA,N), an M by N matrix.
c    On input, the matrix whose decomposition is to be computed.
c    In a least squares data fitting problem, A(I,J) is the
c    value of the J-th basis (model) function at the I-th data point.
c    On output, A contains the output from DQRDC.  The triangular matrix R
c    of the QR factorization is contained in the upper triangle and
c    information needed to recover the orthogonal matrix Q is stored
c    below the diagonal in A and in the vector QRAUX.
c
c    Input, integer LDA, the leading dimension of A.
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, real*8 TOL, a relative tolerance used to determine the
c    numerical rank.  The problem should be scaled so that all the elements
c    of A have roughly the same absolute accuracy EPS.  Then a reasonable
c    value for TOL is roughly EPS divided by the magnitude of the largest
c    element.
c
c    Output, integer KR, the numerical rank.
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 linear
c    system.
c
c    Output, real*8 R(M), the residual, B - A*X.  R may
c    overwrite B.
c
c    Workspace, integer JPVT(N), required if ITASK = 1.
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.  ABS ( A(1,1) ) / ABS ( A(KR,KR) ) is an estimate
c    of the condition number of the matrix of independent columns,
c    and of R.  This estimate will be .le. 1/TOL.
c
c    Workspace, real*8 QRAUX(N), required if ITASK = 1.
c
c    Workspace, real*8 WORK(N), required if ITASK = 1.
c
c    Input, integer ITASK.
c    1, DQRLS factors the matrix A and solves the least squares problem.
c    2, DQRLS assumes that the matrix A was factored with an earlier
c       call to DQRLS, and only solves the least squares problem.
c
c    Output, integer IND, error code.
c    0:  no error
c    -1: LDA .lt. M   (fatal error)
c    -2: N .lt. 1     (fatal error)
c    -3: ITASK .lt. 1 (fatal error)
c
        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 ind
        integer itask
        integer jpvt(n)
        integer kr
        real*8 qraux(n)
        real*8 r(m)
        real*8 tol
        real*8 work(n)
        real*8 x(n)
      
        if ( lda .lt. m ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'DQRLS - Fatal error!'
          write ( *, '(a)' ) '  LDA .lt. M.'
          stop
        end if
      
        if ( n .le. 0 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'DQRLS - Fatal error!'
          write ( *, '(a)' ) '  N .le. 0.'
          stop
        end if
      
        if ( itask .lt. 1 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'DQRLS - Fatal error!'
          write ( *, '(a)' ) '  ITASK .lt. 1.'
          stop
        end if
      
        ind = 0

c  Factor the matrix.

        if ( itask .eq. 1 ) then
          call dqrank ( a, lda, m, n, tol, kr, jpvt, qraux, work )
        end if

c  Solve the least-squares problem.

        call dqrlss ( a, lda, m, n, kr, b, x, r, jpvt, qraux )
      
        return
      end
      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
      subroutine dqrsl ( a, lda, n, k,qraux,y,qy,qty,b,rsd,ab,job,info)
      
c*****************************************************************************80
c
c! DQRSL computes transformations, projections, and least squares solutions.
c
c  Discussion:
c
c    DQRSL requires the output of DQRDC.
c
c    For K .le. min(N,P), let AK be the matrix
c
c      AK = ( A(JPVT(1)), A(JPVT(2)), ..., A(JPVT(K)) )
c
c    formed from columns JPVT(1), ..., JPVT(K) of the original
c    N by P matrix A that was input to DQRDC.  If no pivoting was
c    done, AK consists of the first K columns of A in their
c    original order.  DQRDC produces a factored orthogonal matrix Q
c    and an upper triangular matrix R such that
c
c      AK = Q * (R)
c               (0)
c
c    This information is contained in coded form in the arrays
c    A and QRAUX.
c
c    The parameters QY, QTY, B, RSD, and AB are not referenced
c    if their computation is not requested and in this case
c    can be replaced by dummy variables in the calling program.
c    To save storage, the user may in some cases use the same
c    array for different parameters in the calling sequence.  A
c    frequently occuring example is when one wishes to compute
c    any of B, RSD, or AB and does not need Y or QTY.  In this
c    case one may identify Y, QTY, and one of B, RSD, or AB, while
c    providing separate arrays for anything else that is to be
c    computed.
c
c    Thus the calling sequence
c
c      call dqrsl ( a, lda, n, k, qraux, y, dum, y, b, y, dum, 110, info )
c
c    will result in the computation of B and RSD, with RSD
c    overwriting Y.  More generally, each item in the following
c    list contains groups of permissible identifications for
c    a single calling sequence.
c
c      1. (Y,QTY,B) (RSD) (AB) (QY)
c
c      2. (Y,QTY,RSD) (B) (AB) (QY)
c
c      3. (Y,QTY,AB) (B) (RSD) (QY)
c
c      4. (Y,QY) (QTY,B) (RSD) (AB)
c
c      5. (Y,QY) (QTY,RSD) (B) (AB)
c
c      6. (Y,QY) (QTY,AB) (B) (RSD)
c
c    In any group the value returned in the array allocated to
c    the group corresponds to the last member of the group.
c
c  Modified:
c
c    15 April 2012
c
c  Author:
c
c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart
c
c  Reference:
c
c    Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c  Parameters:
c
c    Input, real*8 A(LDA,P), contains the output of DQRDC.
c
c    Input, integer LDA, the leading dimension of the array A.
c
c    Input, integer N, the number of rows of the matrix AK.  It 
c    must have the same value as N in DQRDC.
c
c    Input, integer K, the number of columns of the matrix AK.  K
c    must not be greater than min(N,P), where P is the same as in the
c    calling sequence to DQRDC.
c
c    Input, real*8 QRAUX(P), the auxiliary output from DQRDC.
c
c    Input, real*8 Y(N), a vector to be manipulated by DQRSL.
c
c    Output, real*8 QY(N), contains Q * Y, if requested.
c
c    Output, real*8 QTY(N), contains Q' * Y, if requested.
c
c    Output, real*8 B(K), the solution of the least squares problem
c      minimize norm2 ( Y - AK * B),
c    if its computation has been requested.  Note that if pivoting was
c    requested in DQRDC, the J-th component of B will be associated with
c    column JPVT(J) of the original matrix A that was input into DQRDC.
c
c    Output, real*8 RSD(N), the least squares residual Y - AK * B,
c    if its computation has been requested.  RSD is also the orthogonal
c    projection of Y onto the orthogonal complement of the column space
c    of AK.
c
c    Output, real*8 AB(N), the least squares approximation Ak * B,
c    if its computation has been requested.  AB is also the orthogonal
c    projection of Y onto the column space of A.
c
c    Input, integer JOB, specifies what is to be computed.  JOB has
c    the decimal expansion ABCDE, with the following meaning:
c
c      if A /= 0, compute QY.
c      if B /= 0, compute QTY.
c      if C /= 0, compute QTY and B.
c      if D /= 0, compute QTY and RSD.
c      if E /= 0, compute QTY and AB.
c
c    Note that a request to compute B, RSD, or AB automatically triggers
c    the computation of QTY, for which an array must be provided in the
c    calling sequence.
c
c    Output, integer INFO, is zero unless the computation of B has
c    been requested and R is exactly singular.  In this case, INFO is the
c    index of the first zero diagonal element of R, and B is left unaltered.

        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer k
        integer lda
        integer n
      
        real*8 a(lda,*)
        real*8 ab(n)
        real*8 b(k)
        logical cab
        logical cb
        logical cqty
        logical cqy
        logical cr
        integer info
        integer j
        integer jj
        integer job
        integer ju
        integer kp1
        real*8 qraux(*)
        real*8 qty(n)
        real*8 qy(n)
        real*8 rsd(n)
        real*8 ddot
        real*8 t
        real*8 temp
        real*8 y(n)

c  set info flag.

        info = 0

c  Determine what is to be computed.

        cqy =        job / 10000         /= 0
        cqty = mod ( job,  10000 )       /= 0
        cb =   mod ( job,   1000 ) / 100 /= 0
        cr =   mod ( job,    100 ) /  10 /= 0
        cab =  mod ( job,     10 )       /= 0
      
        ju = min ( k, n-1 )

c  Special action when N = 1.

        if ( ju .eq. 0 ) then
      
          if ( cqy ) then
            qy(1) = y(1)
          end if
      
          if ( cqty ) then
            qty(1) = y(1)
          end if
      
          if ( cab ) then
            ab(1) = y(1)
           end if
      
          if ( cb ) then
      
            if ( a(1,1) .eq. 0.0D+00 ) then
              info = 1
            else
              b(1) = y(1) / a(1,1)
            end if
      
          end if
      
          if ( cr ) then
            rsd(1) = 0.0D+00
          end if
      
          return
      
        end if

c  Set up to compute QY or QTY.

        if ( cqy ) then
          qy(1:n) = y(1:n)
        end if
      
        if ( cqty ) then
          qty(1:n) = y(1:n)
        end if

c  Compute QY.

        if ( cqy ) then
      
          do jj = 1, ju
      
            j = ju - jj + 1
      
            if ( qraux(j) /= 0.0D+00 ) then
              temp = a(j,j)
              a(j,j) = qraux(j)
              t = - ddot ( n-j+1, a(j,j), 1, qy(j), 1 ) / a(j,j)
              call daxpy ( n-j+1, t, a(j,j), 1, qy(j), 1 )
              a(j,j) = temp
            end if
      
          end do
      
        end if

c  Compute Q'*Y.

           if ( cqty ) then
      
              do j = 1, ju
                 if ( qraux(j) /= 0.0D+00 ) then
                    temp = a(j,j)
                    a(j,j) = qraux(j)
                    t = - ddot ( n-j+1, a(j,j), 1, qty(j), 1 ) / a(j,j)
                    call daxpy ( n-j+1, t, a(j,j), 1, qty(j), 1 )
                    a(j,j) = temp
                 end if
              end do
      
           end if

c  Set up to compute B, RSD, or AB.

           if ( cb ) then
             b(1:k) = qty(1:k)
           end if
      
           kp1 = k + 1
      
           if ( cab ) then
             ab(1:k) = qty(1:k)
           end if
      
           if ( cr .and. (k .lt. n) ) then
             rsd(k+1:n) = qty(k+1:n)
           end if
      
           if ( cab .and. (k+1 .le. n) ) then
              ab(k+1:n) = 0.0D+00
           end if
      
           if ( cr ) then
              rsd(1:k) = 0.0D+00
           end if

c  Compute B.

           if ( cb ) then
      
              do jj = 1, k
      
                 j = k - jj + 1
      
                 if ( a(j,j) .eq. 0.0D+00 ) then
                    info = j
c                    exit
                 go to 10
                 end if
      
                 b(j) = b(j)/a(j,j)
      
                 if ( j /= 1 ) then
                    t = -b(j)
                    call daxpy ( j-1, t, a(1,j), 1, b, 1 )
                 end if
      
              end do
   10      continue   
           end if
      
           if ( cr .or. cab ) then

c  Compute RSD or AB as required.

              do jj = 1, ju
      
                 j = ju - jj + 1
      
                 if ( qraux(j) /= 0.0D+00 ) then
      
                    temp = a(j,j)
                    a(j,j) = qraux(j)
      
                    if ( cr ) then
                     t = - ddot ( n-j+1, a(j,j), 1, rsd(j), 1 ) / a(j,j)
                     call daxpy ( n-j+1, t, a(j,j), 1, rsd(j), 1 )
                    end if
      
                    if ( cab ) then
                      t = - ddot ( n-j+1, a(j,j), 1, ab(j), 1 ) / a(j,j)
                       call daxpy ( n-j+1, t, a(j,j), 1, ab(j), 1 )
                    end if
      
                    a(j,j) = temp
      
                 end if
      
              end do
      
        end if
      
        return
      end
      
      subroutine dscal ( n, sa, x, incx )
      
c*****************************************************************************80
c
c! DSCAL scales a vector by a constant.
c
c  Discussion:
c
c    This routine uses double precision real arithmetic.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    08 April 1999
c
c  Author:
c
c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
c    David Kincaid, Fred Krogh.
c    FORTRAN90 version by John Burkardt.
c
c  Reference:
c
c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
c    Algorithm 539, 
c    Basic Linear Algebra Subprograms for Fortran Usage,
c    ACM Transactions on Mathematical Software,
c    Volume 5, Number 3, September 1979, pages 308-323.
c
c  Parameters:
c
c    Input, integer N, the number of entries in the vector.
c
c    Input, real*8 SA, the multiplier.
c
c    Input/output, real*8 X(*), the vector to be scaled.
c
c    Input, integer INCX, the increment between successive 
c    entries of X.
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer i
        integer incx
        integer ix
        integer m
        integer n
        real*8 sa
        real*8 x(*)
      
        if ( n .le. 0 ) then
      
        else if ( incx .eq. 1 ) then
      
          m = mod ( n, 5 )
      
          x(1:m) = sa * x(1:m)
      
          do i = m+1, n, 5
            x(i)   = sa * x(i)
            x(i+1) = sa * x(i+1)
            x(i+2) = sa * x(i+2)
            x(i+3) = sa * x(i+3)
            x(i+4) = sa * x(i+4)
          end do
      
        else
      
          if ( 0 .le. incx ) then
            ix = 1
          else
            ix = ( - n + 1 ) * incx + 1
          end if
      
          do i = 1, n
            x(ix) = sa * x(ix)
            ix = ix + incx
          end do
      
        end if
      
        return
      end
      
      subroutine dswap ( n, x, incx, y, incy )
      
c*****************************************************************************80
c
c! DSWAP interchanges two vectors.
c
c  Discussion:
c
c    This routine uses double precision real arithmetic.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license. 
c
c  Modified:
c
c    08 April 1999
c
c  Author:
c
c    Original FORTRAN77 version by Charles Lawson, Richard Hanson, 
c    David Kincaid, Fred Krogh.
c    FORTRAN90 version by John Burkardt.
c
c  Reference:
c
c    Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,
c    LINPACK User's Guide,
c    SIAM, 1979,
c    ISBN13: 978-0-898711-72-1,
c    LC: QA214.L56.
c
c    Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,
c    Algorithm 539, 
c    Basic Linear Algebra Subprograms for Fortran Usage,
c    ACM Transactions on Mathematical Software, 
c    Volume 5, Number 3, September 1979, pages 308-323.
c
c  Parameters:
c
c    Input, integer N, the number of entries in the vectors.
c
c    Input/output, real*8 X(*), one of the vectors to swap.
c
c    Input, integer INCX, the increment between successive 
c    entries of X.
c
c    Input/output, real*8 Y(*), one of the vectors to swap.
c
c    Input, integer INCY, the increment between successive 
c    elements of Y.
c
        implicit real*8 (a-h,o-z)
        implicit integer (i-n)
      
        integer i
        integer incx
        integer incy
        integer ix
        integer iy
        integer m
        integer n
        real*8 temp
        real*8 x(*)
        real*8 y(*)
      
        if ( n .le. 0 ) then
      
        else if ( (incx .eq. 1) .and.( incy .eq. 1) ) then
      
          m = mod ( n, 3 )
      
          do i = 1, m
            temp = x(i)
            x(i) = y(i)
            y(i) = temp
          end do
      
          do i = m + 1, n, 3
      
            temp = x(i)
            x(i) = y(i)
            y(i) = temp
      
            temp = x(i+1)
            x(i+1) = y(i+1)
            y(i+1) = temp
      
            temp = x(i+2)
            x(i+2) = y(i+2)
            y(i+2) = temp
      
          end do
      
        else
      
          if ( 0 .le. incx ) then
            ix = 1
          else
            ix = ( - n + 1 ) * incx + 1
          end if
      
          if ( 0 .le. incy ) then
            iy = 1
          else
            iy = ( - n + 1 ) * incy + 1
          end if
      
          do i = 1, n
            temp = x(ix)
            x(ix) = y(iy)
            y(iy) = temp
            ix = ix + incx
            iy = iy + incy
          end do
      
        end if
      
        return
      end