![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZUNGBR( VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO ) 2: * 3: * -- LAPACK routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 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( * ) 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 (MAX(1,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: * If LWORK = -1, then a workspace query is assumed; the routine 88: * only calculates the optimal size of the WORK array, returns 89: * this value as the first entry of the WORK array, and no error 90: * message related to LWORK is issued by XERBLA. 91: * 92: * INFO (output) INTEGER 93: * = 0: successful exit 94: * < 0: if INFO = -i, the i-th argument had an illegal value 95: * 96: * ===================================================================== 97: * 98: * .. Parameters .. 99: COMPLEX*16 ZERO, ONE 100: PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), 101: $ ONE = ( 1.0D+0, 0.0D+0 ) ) 102: * .. 103: * .. Local Scalars .. 104: LOGICAL LQUERY, WANTQ 105: INTEGER I, IINFO, J, LWKOPT, MN, NB 106: * .. 107: * .. External Functions .. 108: LOGICAL LSAME 109: INTEGER ILAENV 110: EXTERNAL LSAME, ILAENV 111: * .. 112: * .. External Subroutines .. 113: EXTERNAL XERBLA, ZUNGLQ, ZUNGQR 114: * .. 115: * .. Intrinsic Functions .. 116: INTRINSIC MAX, MIN 117: * .. 118: * .. Executable Statements .. 119: * 120: * Test the input arguments 121: * 122: INFO = 0 123: WANTQ = LSAME( VECT, 'Q' ) 124: MN = MIN( M, N ) 125: LQUERY = ( LWORK.EQ.-1 ) 126: IF( .NOT.WANTQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN 127: INFO = -1 128: ELSE IF( M.LT.0 ) THEN 129: INFO = -2 130: ELSE IF( N.LT.0 .OR. ( WANTQ .AND. ( N.GT.M .OR. N.LT.MIN( M, 131: $ K ) ) ) .OR. ( .NOT.WANTQ .AND. ( M.GT.N .OR. M.LT. 132: $ MIN( N, K ) ) ) ) THEN 133: INFO = -3 134: ELSE IF( K.LT.0 ) THEN 135: INFO = -4 136: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 137: INFO = -6 138: ELSE IF( LWORK.LT.MAX( 1, MN ) .AND. .NOT.LQUERY ) THEN 139: INFO = -9 140: END IF 141: * 142: IF( INFO.EQ.0 ) THEN 143: IF( WANTQ ) THEN 144: NB = ILAENV( 1, 'ZUNGQR', ' ', M, N, K, -1 ) 145: ELSE 146: NB = ILAENV( 1, 'ZUNGLQ', ' ', M, N, K, -1 ) 147: END IF 148: LWKOPT = MAX( 1, MN )*NB 149: WORK( 1 ) = LWKOPT 150: END IF 151: * 152: IF( INFO.NE.0 ) THEN 153: CALL XERBLA( 'ZUNGBR', -INFO ) 154: RETURN 155: ELSE IF( LQUERY ) THEN 156: RETURN 157: END IF 158: * 159: * Quick return if possible 160: * 161: IF( M.EQ.0 .OR. N.EQ.0 ) THEN 162: WORK( 1 ) = 1 163: RETURN 164: END IF 165: * 166: IF( WANTQ ) THEN 167: * 168: * Form Q, determined by a call to ZGEBRD to reduce an m-by-k 169: * matrix 170: * 171: IF( M.GE.K ) THEN 172: * 173: * If m >= k, assume m >= n >= k 174: * 175: CALL ZUNGQR( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO ) 176: * 177: ELSE 178: * 179: * If m < k, assume m = n 180: * 181: * Shift the vectors which define the elementary reflectors one 182: * column to the right, and set the first row and column of Q 183: * to those of the unit matrix 184: * 185: DO 20 J = M, 2, -1 186: A( 1, J ) = ZERO 187: DO 10 I = J + 1, M 188: A( I, J ) = A( I, J-1 ) 189: 10 CONTINUE 190: 20 CONTINUE 191: A( 1, 1 ) = ONE 192: DO 30 I = 2, M 193: A( I, 1 ) = ZERO 194: 30 CONTINUE 195: IF( M.GT.1 ) THEN 196: * 197: * Form Q(2:m,2:m) 198: * 199: CALL ZUNGQR( M-1, M-1, M-1, A( 2, 2 ), LDA, TAU, WORK, 200: $ LWORK, IINFO ) 201: END IF 202: END IF 203: ELSE 204: * 205: * Form P', determined by a call to ZGEBRD to reduce a k-by-n 206: * matrix 207: * 208: IF( K.LT.N ) THEN 209: * 210: * If k < n, assume k <= m <= n 211: * 212: CALL ZUNGLQ( M, N, K, A, LDA, TAU, WORK, LWORK, IINFO ) 213: * 214: ELSE 215: * 216: * If k >= n, assume m = n 217: * 218: * Shift the vectors which define the elementary reflectors one 219: * row downward, and set the first row and column of P' to 220: * those of the unit matrix 221: * 222: A( 1, 1 ) = ONE 223: DO 40 I = 2, N 224: A( I, 1 ) = ZERO 225: 40 CONTINUE 226: DO 60 J = 2, N 227: DO 50 I = J - 1, 2, -1 228: A( I, J ) = A( I-1, J ) 229: 50 CONTINUE 230: A( 1, J ) = ZERO 231: 60 CONTINUE 232: IF( N.GT.1 ) THEN 233: * 234: * Form P'(2:n,2:n) 235: * 236: CALL ZUNGLQ( N-1, N-1, N-1, A( 2, 2 ), LDA, TAU, WORK, 237: $ LWORK, IINFO ) 238: END IF 239: END IF 240: END IF 241: WORK( 1 ) = LWKOPT 242: RETURN 243: * 244: * End of ZUNGBR 245: * 246: END