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 / ztrevc.f < prev    next >
Text File  |  1996-09-28  |  13KB  |  385 lines
  1.       SUBROUTINE ZTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
  2.      $                   LDVR, MM, M, WORK, RWORK, INFO )
  3. *
  4. *  -- LAPACK routine (version 2.0) --
  5. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  6. *     Courant Institute, Argonne National Lab, and Rice University
  7. *     September 30, 1994
  8. *
  9. *     .. Scalar Arguments ..
  10.       CHARACTER          HOWMNY, SIDE
  11.       INTEGER            INFO, LDT, LDVL, LDVR, M, MM, N
  12. *     ..
  13. *     .. Array Arguments ..
  14.       LOGICAL            SELECT( * )
  15.       DOUBLE PRECISION   RWORK( * )
  16.       COMPLEX*16         T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
  17.      $                   WORK( * )
  18. *     ..
  19. *
  20. *  Purpose
  21. *  =======
  22. *
  23. *  ZTREVC computes some or all of the right and/or left eigenvectors of
  24. *  a complex upper triangular matrix T.
  25. *
  26. *  The right eigenvector x and the left eigenvector y of T corresponding
  27. *  to an eigenvalue w are defined by:
  28. *
  29. *               T*x = w*x,     y'*T = w*y'
  30. *
  31. *  where y' denotes the conjugate transpose of the vector y.
  32. *
  33. *  If all eigenvectors are requested, the routine may either return the
  34. *  matrices X and/or Y of right or left eigenvectors of T, or the
  35. *  products Q*X and/or Q*Y, where Q is an input unitary
  36. *  matrix. If T was obtained from the Schur factorization of an
  37. *  original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of
  38. *  right or left eigenvectors of A.
  39. *
  40. *  Arguments
  41. *  =========
  42. *
  43. *  SIDE    (input) CHARACTER*1
  44. *          = 'R':  compute right eigenvectors only;
  45. *          = 'L':  compute left eigenvectors only;
  46. *          = 'B':  compute both right and left eigenvectors.
  47. *
  48. *  HOWMNY  (input) CHARACTER*1
  49. *          = 'A':  compute all right and/or left eigenvectors;
  50. *          = 'B':  compute all right and/or left eigenvectors,
  51. *                  and backtransform them using the input matrices
  52. *                  supplied in VR and/or VL;
  53. *          = 'S':  compute selected right and/or left eigenvectors,
  54. *                  specified by the logical array SELECT.
  55. *
  56. *  SELECT  (input) LOGICAL array, dimension (N)
  57. *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be
  58. *          computed.
  59. *          If HOWMNY = 'A' or 'B', SELECT is not referenced.
  60. *          To select the eigenvector corresponding to the j-th
  61. *          eigenvalue, SELECT(j) must be set to .TRUE..
  62. *
  63. *  N       (input) INTEGER
  64. *          The order of the matrix T. N >= 0.
  65. *
  66. *  T       (input/output) COMPLEX*16 array, dimension (LDT,N)
  67. *          The upper triangular matrix T.  T is modified, but restored
  68. *          on exit.
  69. *
  70. *  LDT     (input) INTEGER
  71. *          The leading dimension of the array T. LDT >= max(1,N).
  72. *
  73. *  VL      (input/output) COMPLEX*16 array, dimension (LDVL,MM)
  74. *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
  75. *          contain an N-by-N matrix Q (usually the unitary matrix Q of
  76. *          Schur vectors returned by ZHSEQR).
  77. *          On exit, if SIDE = 'L' or 'B', VL contains:
  78. *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
  79. *          if HOWMNY = 'B', the matrix Q*Y;
  80. *          if HOWMNY = 'S', the left eigenvectors of T specified by
  81. *                           SELECT, stored consecutively in the columns
  82. *                           of VL, in the same order as their
  83. *                           eigenvalues.
  84. *          If SIDE = 'R', VL is not referenced.
  85. *
  86. *  LDVL    (input) INTEGER
  87. *          The leading dimension of the array VL.  LDVL >= max(1,N) if
  88. *          SIDE = 'L' or 'B'; LDVL >= 1 otherwise.
  89. *
  90. *  VR      (input/output) COMPLEX*16 array, dimension (LDVR,MM)
  91. *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
  92. *          contain an N-by-N matrix Q (usually the unitary matrix Q of
  93. *          Schur vectors returned by ZHSEQR).
  94. *          On exit, if SIDE = 'R' or 'B', VR contains:
  95. *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;
  96. *          if HOWMNY = 'B', the matrix Q*X;
  97. *          if HOWMNY = 'S', the right eigenvectors of T specified by
  98. *                           SELECT, stored consecutively in the columns
  99. *                           of VR, in the same order as their
  100. *                           eigenvalues.
  101. *          If SIDE = 'L', VR is not referenced.
  102. *
  103. *  LDVR    (input) INTEGER
  104. *          The leading dimension of the array VR.  LDVR >= max(1,N) if
  105. *           SIDE = 'R' or 'B'; LDVR >= 1 otherwise.
  106. *
  107. *  MM      (input) INTEGER
  108. *          The number of columns in the arrays VL and/or VR. MM >= M.
  109. *
  110. *  M       (output) INTEGER
  111. *          The number of columns in the arrays VL and/or VR actually
  112. *          used to store the eigenvectors.  If HOWMNY = 'A' or 'B', M
  113. *          is set to N.  Each selected eigenvector occupies one
  114. *          column.
  115. *
  116. *  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
  117. *
  118. *  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
  119. *
  120. *  INFO    (output) INTEGER
  121. *          = 0:  successful exit
  122. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  123. *
  124. *  Further Details
  125. *  ===============
  126. *
  127. *  The algorithm used in this program is basically backward (forward)
  128. *  substitution, with scaling to make the the code robust against
  129. *  possible overflow.
  130. *
  131. *  Each eigenvector is normalized so that the element of largest
  132. *  magnitude has magnitude 1; here the magnitude of a complex number
  133. *  (x,y) is taken to be |x| + |y|.
  134. *
  135. *  =====================================================================
  136. *
  137. *     .. Parameters ..
  138.       DOUBLE PRECISION   ZERO, ONE
  139.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  140.       COMPLEX*16         CMZERO, CMONE
  141.       PARAMETER          ( CMZERO = ( 0.0D+0, 0.0D+0 ),
  142.      $                   CMONE = ( 1.0D+0, 0.0D+0 ) )
  143. *     ..
  144. *     .. Local Scalars ..
  145.       LOGICAL            ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
  146.       INTEGER            I, II, IS, J, K, KI
  147.       DOUBLE PRECISION   OVFL, REMAX, SCALE, SMIN, SMLNUM, ULP, UNFL
  148.       COMPLEX*16         CDUM
  149. *     ..
  150. *     .. External Functions ..
  151.       LOGICAL            LSAME
  152.       INTEGER            IZAMAX
  153.       DOUBLE PRECISION   DLAMCH, DZASUM
  154.       EXTERNAL           LSAME, IZAMAX, DLAMCH, DZASUM
  155. *     ..
  156. *     .. External Subroutines ..
  157.       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZLATRS
  158. *     ..
  159. *     .. Intrinsic Functions ..
  160.       INTRINSIC          ABS, DBLE, DCMPLX, DCONJG, DIMAG, MAX
  161. *     ..
  162. *     .. Statement Functions ..
  163.       DOUBLE PRECISION   CABS1
  164. *     ..
  165. *     .. Statement Function definitions ..
  166.       CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
  167. *     ..
  168. *     .. Executable Statements ..
  169. *
  170. *     Decode and test the input parameters
  171. *
  172.       BOTHV = LSAME( SIDE, 'B' )
  173.       RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
  174.       LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
  175. *
  176.       ALLV = LSAME( HOWMNY, 'A' )
  177.       OVER = LSAME( HOWMNY, 'B' ) .OR. LSAME( HOWMNY, 'O' )
  178.       SOMEV = LSAME( HOWMNY, 'S' )
  179. *
  180. *     Set M to the number of columns required to store the selected
  181. *     eigenvectors.
  182. *
  183.       IF( SOMEV ) THEN
  184.          M = 0
  185.          DO 10 J = 1, N
  186.             IF( SELECT( J ) )
  187.      $         M = M + 1
  188.    10    CONTINUE
  189.       ELSE
  190.          M = N
  191.       END IF
  192. *
  193.       INFO = 0
  194.       IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
  195.          INFO = -1
  196.       ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
  197.          INFO = -2
  198.       ELSE IF( N.LT.0 ) THEN
  199.          INFO = -4
  200.       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
  201.          INFO = -6
  202.       ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
  203.          INFO = -8
  204.       ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
  205.          INFO = -10
  206.       ELSE IF( MM.LT.M ) THEN
  207.          INFO = -11
  208.       END IF
  209.       IF( INFO.NE.0 ) THEN
  210.          CALL XERBLA( 'ZTREVC', -INFO )
  211.          RETURN
  212.       END IF
  213. *
  214. *     Quick return if possible.
  215. *
  216.       IF( N.EQ.0 )
  217.      $   RETURN
  218. *
  219. *     Set the constants to control overflow.
  220. *
  221.       UNFL = DLAMCH( 'Safe minimum' )
  222.       OVFL = ONE / UNFL
  223.       CALL DLABAD( UNFL, OVFL )
  224.       ULP = DLAMCH( 'Precision' )
  225.       SMLNUM = UNFL*( N / ULP )
  226. *
  227. *     Store the diagonal elements of T in working array WORK.
  228. *
  229.       DO 20 I = 1, N
  230.          WORK( I+N ) = T( I, I )
  231.    20 CONTINUE
  232. *
  233. *     Compute 1-norm of each column of strictly upper triangular
  234. *     part of T to control overflow in triangular solver.
  235. *
  236.       RWORK( 1 ) = ZERO
  237.       DO 30 J = 2, N
  238.          RWORK( J ) = DZASUM( J-1, T( 1, J ), 1 )
  239.    30 CONTINUE
  240. *
  241.       IF( RIGHTV ) THEN
  242. *
  243. *        Compute right eigenvectors.
  244. *
  245.          IS = M
  246.          DO 80 KI = N, 1, -1
  247. *
  248.             IF( SOMEV ) THEN
  249.                IF( .NOT.SELECT( KI ) )
  250.      $            GO TO 80
  251.             END IF
  252.             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
  253. *
  254.             WORK( 1 ) = CMONE
  255. *
  256. *           Form right-hand side.
  257. *
  258.             DO 40 K = 1, KI - 1
  259.                WORK( K ) = -T( K, KI )
  260.    40       CONTINUE
  261. *
  262. *           Solve the triangular system:
  263. *              (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK.
  264. *
  265.             DO 50 K = 1, KI - 1
  266.                T( K, K ) = T( K, K ) - T( KI, KI )
  267.                IF( CABS1( T( K, K ) ).LT.SMIN )
  268.      $            T( K, K ) = SMIN
  269.    50       CONTINUE
  270. *
  271.             IF( KI.GT.1 ) THEN
  272.                CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', 'Y',
  273.      $                      KI-1, T, LDT, WORK( 1 ), SCALE, RWORK,
  274.      $                      INFO )
  275.                WORK( KI ) = SCALE
  276.             END IF
  277. *
  278. *           Copy the vector x or Q*x to VR and normalize.
  279. *
  280.             IF( .NOT.OVER ) THEN
  281.                CALL ZCOPY( KI, WORK( 1 ), 1, VR( 1, IS ), 1 )
  282. *
  283.                II = IZAMAX( KI, VR( 1, IS ), 1 )
  284.                REMAX = ONE / CABS1( VR( II, IS ) )
  285.                CALL ZDSCAL( KI, REMAX, VR( 1, IS ), 1 )
  286. *
  287.                DO 60 K = KI + 1, N
  288.                   VR( K, IS ) = CMZERO
  289.    60          CONTINUE
  290.             ELSE
  291.                IF( KI.GT.1 )
  292.      $            CALL ZGEMV( 'N', N, KI-1, CMONE, VR, LDVR, WORK( 1 ),
  293.      $                        1, DCMPLX( SCALE ), VR( 1, KI ), 1 )
  294. *
  295.                II = IZAMAX( N, VR( 1, KI ), 1 )
  296.                REMAX = ONE / CABS1( VR( II, KI ) )
  297.                CALL ZDSCAL( N, REMAX, VR( 1, KI ), 1 )
  298.             END IF
  299. *
  300. *           Set back the original diagonal elements of T.
  301. *
  302.             DO 70 K = 1, KI - 1
  303.                T( K, K ) = WORK( K+N )
  304.    70       CONTINUE
  305. *
  306.             IS = IS - 1
  307.    80    CONTINUE
  308.       END IF
  309. *
  310.       IF( LEFTV ) THEN
  311. *
  312. *        Compute left eigenvectors.
  313. *
  314.          IS = 1
  315.          DO 130 KI = 1, N
  316. *
  317.             IF( SOMEV ) THEN
  318.                IF( .NOT.SELECT( KI ) )
  319.      $            GO TO 130
  320.             END IF
  321.             SMIN = MAX( ULP*( CABS1( T( KI, KI ) ) ), SMLNUM )
  322. *
  323.             WORK( N ) = CMONE
  324. *
  325. *           Form right-hand side.
  326. *
  327.             DO 90 K = KI + 1, N
  328.                WORK( K ) = -DCONJG( T( KI, K ) )
  329.    90       CONTINUE
  330. *
  331. *           Solve the triangular system:
  332. *              (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK.
  333. *
  334.             DO 100 K = KI + 1, N
  335.                T( K, K ) = T( K, K ) - T( KI, KI )
  336.                IF( CABS1( T( K, K ) ).LT.SMIN )
  337.      $            T( K, K ) = SMIN
  338.   100       CONTINUE
  339. *
  340.             IF( KI.LT.N ) THEN
  341.                CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
  342.      $                      'Y', N-KI, T( KI+1, KI+1 ), LDT,
  343.      $                      WORK( KI+1 ), SCALE, RWORK, INFO )
  344.                WORK( KI ) = SCALE
  345.             END IF
  346. *
  347. *           Copy the vector x or Q*x to VL and normalize.
  348. *
  349.             IF( .NOT.OVER ) THEN
  350.                CALL ZCOPY( N-KI+1, WORK( KI ), 1, VL( KI, IS ), 1 )
  351. *
  352.                II = IZAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
  353.                REMAX = ONE / CABS1( VL( II, IS ) )
  354.                CALL ZDSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
  355. *
  356.                DO 110 K = 1, KI - 1
  357.                   VL( K, IS ) = CMZERO
  358.   110          CONTINUE
  359.             ELSE
  360.                IF( KI.LT.N )
  361.      $            CALL ZGEMV( 'N', N, N-KI, CMONE, VL( 1, KI+1 ), LDVL,
  362.      $                        WORK( KI+1 ), 1, DCMPLX( SCALE ),
  363.      $                        VL( 1, KI ), 1 )
  364. *
  365.                II = IZAMAX( N, VL( 1, KI ), 1 )
  366.                REMAX = ONE / CABS1( VL( II, KI ) )
  367.                CALL ZDSCAL( N, REMAX, VL( 1, KI ), 1 )
  368.             END IF
  369. *
  370. *           Set back the original diagonal elements of T.
  371. *
  372.             DO 120 K = KI + 1, N
  373.                T( K, K ) = WORK( K+N )
  374.   120       CONTINUE
  375. *
  376.             IS = IS + 1
  377.   130    CONTINUE
  378.       END IF
  379. *
  380.       RETURN
  381. *
  382. *     End of ZTREVC
  383. *
  384.       END
  385.