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 / zungbr.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  225 lines
  1.       SUBROUTINE ZUNGBR( VECT, M, N, K, A, LDA, TAU, WORK, LWORK, 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. *     September 30, 1994
  7. *
  8. *     .. Scalar Arguments ..
  9.       CHARACTER          VECT
  10.       INTEGER            INFO, K, LDA, LWORK, M, N
  11. *     ..
  12. *     .. Array Arguments ..
  13.       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( LWORK )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  ZUNGBR generates one of the complex unitary matrices Q or P**H
  20. *  determined by ZGEBRD when reducing a complex matrix A to bidiagonal
  21. *  form: A = Q * B * P**H.  Q and P**H are defined as products of
  22. *  elementary reflectors H(i) or G(i) respectively.
  23. *
  24. *  If VECT = 'Q', A is assumed to have been an M-by-K matrix, and Q
  25. *  is of order M:
  26. *  if m >= k, Q = H(1) H(2) . . . H(k) and ZUNGBR returns the first n
  27. *  columns of Q, where m >= n >= k;
  28. *  if m < k, Q = H(1) H(2) . . . H(m-1) and ZUNGBR returns Q as an
  29. *  M-by-M matrix.
  30. *
  31. *  If VECT = 'P', A is assumed to have been a K-by-N matrix, and P**H
  32. *  is of order N:
  33. *  if k < n, P**H = G(k) . . . G(2) G(1) and ZUNGBR returns the first m
  34. *  rows of P**H, where n >= m >= k;
  35. *  if k >= n, P**H = G(n-1) . . . G(2) G(1) and ZUNGBR returns P**H as
  36. *  an N-by-N matrix.
  37. *
  38. *  Arguments
  39. *  =========
  40. *
  41. *  VECT    (input) CHARACTER*1
  42. *          Specifies whether the matrix Q or the matrix P**H is
  43. *          required, as defined in the transformation applied by ZGEBRD:
  44. *          = 'Q':  generate Q;
  45. *          = 'P':  generate P**H.
  46. *
  47. *  M       (input) INTEGER
  48. *          The number of rows of the matrix Q or P**H to be returned.
  49. *          M >= 0.
  50. *
  51. *  N       (input) INTEGER
  52. *          The number of columns of the matrix Q or P**H to be returned.
  53. *          N >= 0.
  54. *          If VECT = 'Q', M >= N >= min(M,K);
  55. *          if VECT = 'P', N >= M >= min(N,K).
  56. *
  57. *  K       (input) INTEGER
  58. *          If VECT = 'Q', the number of columns in the original M-by-K
  59. *          matrix reduced by ZGEBRD.
  60. *          If VECT = 'P', the number of rows in the original K-by-N
  61. *          matrix reduced by ZGEBRD.
  62. *          K >= 0.
  63. *
  64. *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
  65. *          On entry, the vectors which define the elementary reflectors,
  66. *          as returned by ZGEBRD.
  67. *          On exit, the M-by-N matrix Q or P**H.
  68. *
  69. *  LDA     (input) INTEGER
  70. *          The leading dimension of the array A. LDA >= M.
  71. *
  72. *  TAU     (input) COMPLEX*16 array, dimension
  73. *                                (min(M,K)) if VECT = 'Q'
  74. *                                (min(N,K)) if VECT = 'P'
  75. *          TAU(i) must contain the scalar factor of the elementary
  76. *          reflector H(i) or G(i), which determines Q or P**H, as
  77. *          returned by ZGEBRD in its array argument TAUQ or TAUP.
  78. *
  79. *  WORK    (workspace/output) COMPLEX*16 array, dimension (LWORK)
  80. *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  81. *
  82. *  LWORK   (input) INTEGER
  83. *          The dimension of the array WORK. LWORK >= max(1,min(M,N)).
  84. *          For optimum performance LWORK >= min(M,N)*NB, where NB
  85. *          is the optimal blocksize.
  86. *
  87. *  INFO    (output) INTEGER
  88. *          = 0:  successful exit
  89. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  90. *
  91. *  =====================================================================
  92. *
  93. *     .. Parameters ..
  94.       COMPLEX*16         ZERO, ONE
  95.       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ),
  96.      $                   ONE = ( 1.0D+0, 0.0D+0 ) )
  97. *     ..
  98. *     .. Local Scalars ..
  99.       LOGICAL            WANTQ
  100.       INTEGER            I, IINFO, J
  101. *     ..
  102. *     .. External Functions ..
  103.       LOGICAL            LSAME
  104.       EXTERNAL           LSAME
  105. *     ..
  106. *     .. External Subroutines ..
  107.       EXTERNAL           XERBLA, ZUNGLQ, ZUNGQR
  108. *     ..
  109. *     .. Intrinsic Functions ..
  110.       INTRINSIC          MAX, MIN
  111. *     ..
  112. *     .. Executable Statements ..
  113. *
  114. *     Test the input arguments
  115. *
  116.       INFO = 0
  117.       WANTQ = LSAME( VECT, 'Q' )
  118.       IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
  119.          INFO = -1
  120.       ELSE IF( M.LT.0 ) THEN
  121.          INFO = -2
  122.       ELSE IF( N.LT.0 .OR. ( WANTQ .AND. ( N.GT.M .OR. N.LT.MIN( M,
  123.      $         K ) ) ) .OR. ( .NOT.WANTQ .AND. ( M.GT.N .OR. M.LT.
  124.      $         MIN( N, K ) ) ) ) THEN
  125.          INFO = -3
  126.       ELSE IF( K.LT.0 ) THEN
  127.          INFO = -4
  128.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  129.          INFO = -6
  130.       ELSE IF( LWORK.LT.MAX( 1, MIN( M, N ) ) ) THEN
  131.          INFO = -9
  132.       END IF
  133.       IF( INFO.NE.0 ) THEN
  134.          CALL XERBLA( 'ZUNGBR', -INFO )
  135.          RETURN
  136.       END IF
  137. *
  138. *     Quick return if possible
  139. *
  140.       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  141.          WORK( 1 ) = 1
  142.          RETURN
  143.       END IF
  144. *
  145.       IF( WANTQ ) THEN
  146. *
  147. *        Form Q, determined by a call to ZGEBRD to reduce an m-by-k
  148. *        matrix
  149. *
  150.          IF( M.GE.K ) THEN
  151. *
  152. *           If m >= k, assume m >= n >= k
  153. *
  154.             CALL ZUNGQR( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  155. *
  156.          ELSE
  157. *
  158. *           If m < k, assume m = n
  159. *
  160. *           Shift the vectors which define the elementary reflectors one
  161. *           column to the right, and set the first row and column of Q
  162. *           to those of the unit matrix
  163. *
  164.             DO 20 J = M, 2, -1
  165.                A( 1, J ) = ZERO
  166.                DO 10 I = J + 1, M
  167.                   A( I, J ) = A( I, J-1 )
  168.    10          CONTINUE
  169.    20       CONTINUE
  170.             A( 1, 1 ) = ONE
  171.             DO 30 I = 2, M
  172.                A( I, 1 ) = ZERO
  173.    30       CONTINUE
  174.             IF( M.GT.1 ) THEN
  175. *
  176. *              Form Q(2:m,2:m)
  177. *
  178.                CALL ZUNGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK,
  179.      $                      LWORK, IINFO )
  180.             END IF
  181.          END IF
  182.       ELSE
  183. *
  184. *        Form P', determined by a call to ZGEBRD to reduce a k-by-n
  185. *        matrix
  186. *
  187.          IF( K.LT.N ) THEN
  188. *
  189. *           If k < n, assume k <= m <= n
  190. *
  191.             CALL ZUNGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO )
  192. *
  193.          ELSE
  194. *
  195. *           If k >= n, assume m = n
  196. *
  197. *           Shift the vectors which define the elementary reflectors one
  198. *           row downward, and set the first row and column of P' to
  199. *           those of the unit matrix
  200. *
  201.             A( 1, 1 ) = ONE
  202.             DO 40 I = 2, N
  203.                A( I, 1 ) = ZERO
  204.    40       CONTINUE
  205.             DO 60 J = 2, N
  206.                DO 50 I = J - 1, 2, -1
  207.                   A( I, J ) = A( I-1, J )
  208.    50          CONTINUE
  209.                A( 1, J ) = ZERO
  210.    60       CONTINUE
  211.             IF( N.GT.1 ) THEN
  212. *
  213. *              Form P'(2:n,2:n)
  214. *
  215.                CALL ZUNGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK,
  216.      $                      LWORK, IINFO )
  217.             END IF
  218.          END IF
  219.       END IF
  220.       RETURN
  221. *
  222. *     End of ZUNGBR
  223. *
  224.       END
  225.