Sources Esope | Cast3M

dqrsl
C DQRSL     SOURCE    PV090527  23/02/07    11:58:38     11592          
      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