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 / dlarfg.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  139 lines
  1.       SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
  2. *
  3. *  -- LAPACK auxiliary routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INCX, N
  10.       DOUBLE PRECISION   ALPHA, TAU
  11. *     ..
  12. *     .. Array Arguments ..
  13.       DOUBLE PRECISION   X( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DLARFG generates a real elementary reflector H of order n, such
  20. *  that
  21. *
  22. *        H * ( alpha ) = ( beta ),   H' * H = I.
  23. *            (   x   )   (   0  )
  24. *
  25. *  where alpha and beta are scalars, and x is an (n-1)-element real
  26. *  vector. H is represented in the form
  27. *
  28. *        H = I - tau * ( 1 ) * ( 1 v' ) ,
  29. *                      ( v )
  30. *
  31. *  where tau is a real scalar and v is a real (n-1)-element
  32. *  vector.
  33. *
  34. *  If the elements of x are all zero, then tau = 0 and H is taken to be
  35. *  the unit matrix.
  36. *
  37. *  Otherwise  1 <= tau <= 2.
  38. *
  39. *  Arguments
  40. *  =========
  41. *
  42. *  N       (input) INTEGER
  43. *          The order of the elementary reflector.
  44. *
  45. *  ALPHA   (input/output) DOUBLE PRECISION
  46. *          On entry, the value alpha.
  47. *          On exit, it is overwritten with the value beta.
  48. *
  49. *  X       (input/output) DOUBLE PRECISION array, dimension
  50. *                         (1+(N-2)*abs(INCX))
  51. *          On entry, the vector x.
  52. *          On exit, it is overwritten with the vector v.
  53. *
  54. *  INCX    (input) INTEGER
  55. *          The increment between elements of X. INCX > 0.
  56. *
  57. *  TAU     (output) DOUBLE PRECISION
  58. *          The value tau.
  59. *
  60. *  =====================================================================
  61. *
  62. *     .. Parameters ..
  63.       DOUBLE PRECISION   ONE, ZERO
  64.       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  65. *     ..
  66. *     .. Local Scalars ..
  67.       INTEGER            J, KNT
  68.       DOUBLE PRECISION   BETA, RSAFMN, SAFMIN, XNORM
  69. *     ..
  70. *     .. External Functions ..
  71.       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
  72.       EXTERNAL           DLAMCH, DLAPY2, DNRM2
  73. *     ..
  74. *     .. Intrinsic Functions ..
  75.       INTRINSIC          ABS, SIGN
  76. *     ..
  77. *     .. External Subroutines ..
  78.       EXTERNAL           DSCAL
  79. *     ..
  80. *     .. Executable Statements ..
  81. *
  82.       IF( N.LE.1 ) THEN
  83.          TAU = ZERO
  84.          RETURN
  85.       END IF
  86. *
  87.       XNORM = DNRM2( N-1, X, INCX )
  88. *
  89.       IF( XNORM.EQ.ZERO ) THEN
  90. *
  91. *        H  =  I
  92. *
  93.          TAU = ZERO
  94.       ELSE
  95. *
  96. *        general case
  97. *
  98.          BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  99.          SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
  100.          IF( ABS( BETA ).LT.SAFMIN ) THEN
  101. *
  102. *           XNORM, BETA may be inaccurate; scale X and recompute them
  103. *
  104.             RSAFMN = ONE / SAFMIN
  105.             KNT = 0
  106.    10       CONTINUE
  107.             KNT = KNT + 1
  108.             CALL DSCAL( N-1, RSAFMN, X, INCX )
  109.             BETA = BETA*RSAFMN
  110.             ALPHA = ALPHA*RSAFMN
  111.             IF( ABS( BETA ).LT.SAFMIN )
  112.      $         GO TO 10
  113. *
  114. *           New BETA is at most 1, at least SAFMIN
  115. *
  116.             XNORM = DNRM2( N-1, X, INCX )
  117.             BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  118.             TAU = ( BETA-ALPHA ) / BETA
  119.             CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
  120. *
  121. *           If ALPHA is subnormal, it may lose relative accuracy
  122. *
  123.             ALPHA = BETA
  124.             DO 20 J = 1, KNT
  125.                ALPHA = ALPHA*SAFMIN
  126.    20       CONTINUE
  127.          ELSE
  128.             TAU = ( BETA-ALPHA ) / BETA
  129.             CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
  130.             ALPHA = BETA
  131.          END IF
  132.       END IF
  133. *
  134.       RETURN
  135. *
  136. *     End of DLARFG
  137. *
  138.       END
  139.