share/odepack/fortran/dhels.f

   1 *DECK DHELS
   2       SUBROUTINE DHELS (A, LDA, N, Q, B)
   3       INTEGER LDA, N
   4       DOUBLE PRECISION A(LDA,*), B(*), Q(*)
   5 C-----------------------------------------------------------------------
   6 C This is part of the LINPACK routine DGESL with changes
   7 C due to the fact that A is an upper Hessenberg matrix.
   8 C-----------------------------------------------------------------------
   9 C     DHELS solves the least squares problem
  10 C
  11 C           min (b-A*x, b-A*x)
  12 C
  13 C     using the factors computed by DHEQR.
  14 C
  15 C     On entry
  16 C
  17 C        A       DOUBLE PRECISION(LDA, N)
  18 C                the output from DHEQR which contains the upper
  19 C                triangular factor R in the QR decomposition of A.
  20 C
  21 C        LDA     INTEGER
  22 C                the leading dimension of the array  A .
  23 C
  24 C        N       INTEGER
  25 C                A is originally an (N+1) by N matrix.
  26 C
  27 C        Q       DOUBLE PRECISION(2*N)
  28 C                The coefficients of the N givens rotations
  29 C                used in the QR factorization of A.
  30 C
  31 C        B       DOUBLE PRECISION(N+1)
  32 C                the right hand side vector.
  33 C
  34 C     On return
  35 C
  36 C        B       the solution vector  x .
  37 C
  38 C     Modification of LINPACK, by Peter Brown, LLNL.
  39 C     Written 1/13/86.  This version dated 6/20/01.
  40 C
  41 C     BLAS called: DAXPY
  42 C-----------------------------------------------------------------------
  43       INTEGER IQ, K, KB, KP1
  44       DOUBLE PRECISION C, S, T, T1, T2
  45 C
  46 C        Minimize (b-A*x, b-A*x)
  47 C        First form Q*b.
  48 C
  49          DO 20 K = 1, N
  50             KP1 = K + 1
  51             IQ = 2*(K-1) + 1
  52             C = Q(IQ)
  53             S = Q(IQ+1)
  54             T1 = B(K)
  55             T2 = B(KP1)
  56             B(K) = C*T1 - S*T2
  57             B(KP1) = S*T1 + C*T2
  58    20    CONTINUE
  59 C
  60 C        Now solve  R*x = Q*b.
  61 C
  62          DO 40 KB = 1, N
  63             K = N + 1 - KB
  64             B(K) = B(K)/A(K,K)
  65             T = -B(K)
  66             CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)
  67    40    CONTINUE
  68       RETURN
  69 C----------------------- End of Subroutine DHELS -----------------------
  70       END