share/odepack/fortran/dhesl.f

   1 *DECK DHESL
   2       SUBROUTINE DHESL (A, LDA, N, IPVT, B)
   3       INTEGER LDA, N, IPVT(*)
   4       DOUBLE PRECISION A(LDA,*), B(*)
   5 C-----------------------------------------------------------------------
   6 C This is essentially the LINPACK routine DGESL except for changes
   7 C due to the fact that A is an upper Hessenberg matrix.
   8 C-----------------------------------------------------------------------
   9 C     DHESL solves the real system A * x = b
  10 C     using the factors computed by DHEFA.
  11 C
  12 C     On entry
  13 C
  14 C        A       DOUBLE PRECISION(LDA, N)
  15 C                the output from DHEFA.
  16 C
  17 C        LDA     INTEGER
  18 C                the leading dimension of the array  A .
  19 C
  20 C        N       INTEGER
  21 C                the order of the matrix  A .
  22 C
  23 C        IPVT    INTEGER(N)
  24 C                the pivot vector from DHEFA.
  25 C
  26 C        B       DOUBLE PRECISION(N)
  27 C                the right hand side vector.
  28 C
  29 C     On return
  30 C
  31 C        B       the solution vector  x .
  32 C
  33 C     Modification of LINPACK, by Peter Brown, LLNL.
  34 C     Written 7/20/83.  This version dated 6/20/01.
  35 C
  36 C     BLAS called: DAXPY
  37 C-----------------------------------------------------------------------
  38       INTEGER K, KB, L, NM1
  39       DOUBLE PRECISION T
  40 C
  41       NM1 = N - 1
  42 C
  43 C        Solve  A * x = b
  44 C        First solve  L*y = b
  45 C
  46          IF (NM1 .LT. 1) GO TO 30
  47          DO 20 K = 1, NM1
  48             L = IPVT(K)
  49             T = B(L)
  50             IF (L .EQ. K) GO TO 10
  51                B(L) = B(K)
  52                B(K) = T
  53    10       CONTINUE
  54             B(K+1) = B(K+1) + T*A(K+1,K)
  55    20    CONTINUE
  56    30    CONTINUE
  57 C
  58 C        Now solve  U*x = y
  59 C
  60          DO 40 KB = 1, N
  61             K = N + 1 - KB
  62             B(K) = B(K)/A(K,K)
  63             T = -B(K)
  64             CALL DAXPY (K-1, T, A(1,K), 1, B(1), 1)
  65    40    CONTINUE
  66       RETURN
  67 C----------------------- End of Subroutine DHESL -----------------------
  68       END