home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / dgehd2.f < prev    next >
Text File  |  1996-09-28  |  5KB  |  151 lines
  1.       SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO )
  2. *
  3. *  -- LAPACK routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     October 31, 1992
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            IHI, ILO, INFO, LDA, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
  13. *     ..
  14. *
  15. *  Purpose
  16. *  =======
  17. *
  18. *  DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
  19. *  an orthogonal similarity transformation:  Q' * A * Q = H .
  20. *
  21. *  Arguments
  22. *  =========
  23. *
  24. *  N       (input) INTEGER
  25. *          The order of the matrix A.  N >= 0.
  26. *
  27. *  ILO     (input) INTEGER
  28. *  IHI     (input) INTEGER
  29. *          It is assumed that A is already upper triangular in rows
  30. *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
  31. *          set by a previous call to DGEBAL; otherwise they should be
  32. *          set to 1 and N respectively. See Further Details.
  33. *          1 <= ILO <= IHI <= max(1,N).
  34. *
  35. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  36. *          On entry, the n by n general matrix to be reduced.
  37. *          On exit, the upper triangle and the first subdiagonal of A
  38. *          are overwritten with the upper Hessenberg matrix H, and the
  39. *          elements below the first subdiagonal, with the array TAU,
  40. *          represent the orthogonal matrix Q as a product of elementary
  41. *          reflectors. See Further Details.
  42. *
  43. *  LDA     (input) INTEGER
  44. *          The leading dimension of the array A.  LDA >= max(1,N).
  45. *
  46. *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
  47. *          The scalar factors of the elementary reflectors (see Further
  48. *          Details).
  49. *
  50. *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
  51. *
  52. *  INFO    (output) INTEGER
  53. *          = 0:  successful exit.
  54. *          < 0:  if INFO = -i, the i-th argument had an illegal value.
  55. *
  56. *  Further Details
  57. *  ===============
  58. *
  59. *  The matrix Q is represented as a product of (ihi-ilo) elementary
  60. *  reflectors
  61. *
  62. *     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
  63. *
  64. *  Each H(i) has the form
  65. *
  66. *     H(i) = I - tau * v * v'
  67. *
  68. *  where tau is a real scalar, and v is a real vector with
  69. *  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
  70. *  exit in A(i+2:ihi,i), and tau in TAU(i).
  71. *
  72. *  The contents of A are illustrated by the following example, with
  73. *  n = 7, ilo = 2 and ihi = 6:
  74. *
  75. *  on entry,                        on exit,
  76. *
  77. *  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
  78. *  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
  79. *  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
  80. *  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
  81. *  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
  82. *  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
  83. *  (                         a )    (                          a )
  84. *
  85. *  where a denotes an element of the original matrix A, h denotes a
  86. *  modified element of the upper Hessenberg matrix H, and vi denotes an
  87. *  element of the vector defining H(i).
  88. *
  89. *  =====================================================================
  90. *
  91. *     .. Parameters ..
  92.       DOUBLE PRECISION   ONE
  93.       PARAMETER          ( ONE = 1.0D+0 )
  94. *     ..
  95. *     .. Local Scalars ..
  96.       INTEGER            I
  97.       DOUBLE PRECISION   AII
  98. *     ..
  99. *     .. External Subroutines ..
  100.       EXTERNAL           DLARF, DLARFG, XERBLA
  101. *     ..
  102. *     .. Intrinsic Functions ..
  103.       INTRINSIC          MAX, MIN
  104. *     ..
  105. *     .. Executable Statements ..
  106. *
  107. *     Test the input parameters
  108. *
  109.       INFO = 0
  110.       IF( N.LT.0 ) THEN
  111.          INFO = -1
  112.       ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
  113.          INFO = -2
  114.       ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
  115.          INFO = -3
  116.       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  117.          INFO = -5
  118.       END IF
  119.       IF( INFO.NE.0 ) THEN
  120.          CALL XERBLA( 'DGEHD2', -INFO )
  121.          RETURN
  122.       END IF
  123. *
  124.       DO 10 I = ILO, IHI - 1
  125. *
  126. *        Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
  127. *
  128.          CALL DLARFG( IHI-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1,
  129.      $                TAU( I ) )
  130.          AII = A( I+1, I )
  131.          A( I+1, I ) = ONE
  132. *
  133. *        Apply H(i) to A(1:ihi,i+1:ihi) from the right
  134. *
  135.          CALL DLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ),
  136.      $               A( 1, I+1 ), LDA, WORK )
  137. *
  138. *        Apply H(i) to A(i+1:ihi,i+1:n) from the left
  139. *
  140.          CALL DLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1, TAU( I ),
  141.      $               A( I+1, I+1 ), LDA, WORK )
  142. *
  143.          A( I+1, I ) = AII
  144.    10 CONTINUE
  145. *
  146.       RETURN
  147. *
  148. *     End of DGEHD2
  149. *
  150.       END
  151.