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 / zlarfg.f < prev    next >
Text File  |  1996-09-28  |  4KB  |  147 lines
  1.       SUBROUTINE ZLARFG( 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.       COMPLEX*16         ALPHA, TAU
  11. *     ..
  12. *     .. Array Arguments ..
  13.       COMPLEX*16         X( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  ZLARFG generates a complex 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, with beta real, and x is an
  26. *  (n-1)-element complex vector. H is represented in the form
  27. *
  28. *        H = I - tau * ( 1 ) * ( 1 v' ) ,
  29. *                      ( v )
  30. *
  31. *  where tau is a complex scalar and v is a complex (n-1)-element
  32. *  vector. Note that H is not hermitian.
  33. *
  34. *  If the elements of x are all zero and alpha is real, then tau = 0
  35. *  and H is taken to be the unit matrix.
  36. *
  37. *  Otherwise  1 <= real(tau) <= 2  and  abs(tau-1) <= 1 .
  38. *
  39. *  Arguments
  40. *  =========
  41. *
  42. *  N       (input) INTEGER
  43. *          The order of the elementary reflector.
  44. *
  45. *  ALPHA   (input/output) COMPLEX*16
  46. *          On entry, the value alpha.
  47. *          On exit, it is overwritten with the value beta.
  48. *
  49. *  X       (input/output) COMPLEX*16 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) COMPLEX*16
  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   ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM
  69. *     ..
  70. *     .. External Functions ..
  71.       DOUBLE PRECISION   DLAMCH, DLAPY3, DZNRM2
  72.       COMPLEX*16         ZLADIV
  73.       EXTERNAL           DLAMCH, DLAPY3, DZNRM2, ZLADIV
  74. *     ..
  75. *     .. Intrinsic Functions ..
  76.       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, SIGN
  77. *     ..
  78. *     .. External Subroutines ..
  79.       EXTERNAL           ZDSCAL, ZSCAL
  80. *     ..
  81. *     .. Executable Statements ..
  82. *
  83.       IF( N.LE.0 ) THEN
  84.          TAU = ZERO
  85.          RETURN
  86.       END IF
  87. *
  88.       XNORM = DZNRM2( N-1, X, INCX )
  89.       ALPHR = DBLE( ALPHA )
  90.       ALPHI = DIMAG( ALPHA )
  91. *
  92.       IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN
  93. *
  94. *        H  =  I
  95. *
  96.          TAU = ZERO
  97.       ELSE
  98. *
  99. *        general case
  100. *
  101.          BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  102.          SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
  103.          RSAFMN = ONE / SAFMIN
  104. *
  105.          IF( ABS( BETA ).LT.SAFMIN ) THEN
  106. *
  107. *           XNORM, BETA may be inaccurate; scale X and recompute them
  108. *
  109.             KNT = 0
  110.    10       CONTINUE
  111.             KNT = KNT + 1
  112.             CALL ZDSCAL( N-1, RSAFMN, X, INCX )
  113.             BETA = BETA*RSAFMN
  114.             ALPHI = ALPHI*RSAFMN
  115.             ALPHR = ALPHR*RSAFMN
  116.             IF( ABS( BETA ).LT.SAFMIN )
  117.      $         GO TO 10
  118. *
  119. *           New BETA is at most 1, at least SAFMIN
  120. *
  121.             XNORM = DZNRM2( N-1, X, INCX )
  122.             ALPHA = DCMPLX( ALPHR, ALPHI )
  123.             BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  124.             TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
  125.             ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
  126.             CALL ZSCAL( N-1, ALPHA, X, INCX )
  127. *
  128. *           If ALPHA is subnormal, it may lose relative accuracy
  129. *
  130.             ALPHA = BETA
  131.             DO 20 J = 1, KNT
  132.                ALPHA = ALPHA*SAFMIN
  133.    20       CONTINUE
  134.          ELSE
  135.             TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA )
  136.             ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA )
  137.             CALL ZSCAL( N-1, ALPHA, X, INCX )
  138.             ALPHA = BETA
  139.          END IF
  140.       END IF
  141. *
  142.       RETURN
  143. *
  144. *     End of ZLARFG
  145. *
  146.       END
  147.