dhels.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number 
C W-7405-Eng-48.
C
      SUBROUTINE dhels (A, LDA, N, Q, B)
C
C***BEGIN PROLOGUE  DHELS
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C This is similar to the LINPACK routine DGESL except that
C A is an upper Hessenberg matrix.
C
C     DHELS solves the least squares problem
C
C           MIN (B-A*X,B-A*X)
C
C     using the factors computed by DHEQR.
C
C     On entry
C
C        A       DOUBLE PRECISION (LDA, N)
C                The output from DHEQR which contains the upper
C                triangular factor R in the QR decomposition of A.
C
C        LDA     INTEGER
C                The leading dimension of the array  A .
C
C        N       INTEGER
C                A is originally an (N+1) by N matrix.
C
C        Q       DOUBLE PRECISION(2*N)
C                The coefficients of the N givens rotations
C                used in the QR factorization of A.
C
C        B       DOUBLE PRECISION(N+1)
C                The right hand side vector.
C
C
C     On return
C
C        B       The solution vector X.
C
C
C     Modification of LINPACK.
C     Peter Brown, Lawrence Livermore Natl. Lab.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   DAXPY 
C
C***END PROLOGUE  DHELS
C
      INTEGER lda, n
      DOUBLE PRECISION a(lda,*), b(*), q(*)
      INTEGER iq, k, kb, kp1
      DOUBLE PRECISION c, s, t, t1, t2
C
C        Minimize (B-A*X,B-A*X).
C        First form Q*B.
C
         DO 20 k = 1, n
            kp1 = k + 1
            iq = 2*(k-1) + 1
            c = q(iq)
            s = q(iq+1)
            t1 = b(k)
            t2 = b(kp1)
            b(k) = c*t1 - s*t2
            b(kp1) = s*t1 + c*t2
   20    CONTINUE
C
C        Now solve R*X = Q*B.
C
         DO 40 kb = 1, n
            k = n + 1 - kb
            b(k) = b(k)/a(k,k)
            t = -b(k)
            CALL daxpy(k-1, t, a(1,k), 1, b(1), 1)
   40    CONTINUE
      RETURN
C
C------END OF SUBROUTINE DHELS------------------------------------------
      END