![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE ZGGRQF( M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, 2: $ LWORK, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: INTEGER INFO, LDA, LDB, LWORK, M, N, P 11: * .. 12: * .. Array Arguments .. 13: COMPLEX*16 A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ), 14: $ WORK( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * ZGGRQF computes a generalized RQ factorization of an M-by-N matrix A 21: * and a P-by-N matrix B: 22: * 23: * A = R*Q, B = Z*T*Q, 24: * 25: * where Q is an N-by-N unitary matrix, Z is a P-by-P unitary 26: * matrix, and R and T assume one of the forms: 27: * 28: * if M <= N, R = ( 0 R12 ) M, or if M > N, R = ( R11 ) M-N, 29: * N-M M ( R21 ) N 30: * N 31: * 32: * where R12 or R21 is upper triangular, and 33: * 34: * if P >= N, T = ( T11 ) N , or if P < N, T = ( T11 T12 ) P, 35: * ( 0 ) P-N P N-P 36: * N 37: * 38: * where T11 is upper triangular. 39: * 40: * In particular, if B is square and nonsingular, the GRQ factorization 41: * of A and B implicitly gives the RQ factorization of A*inv(B): 42: * 43: * A*inv(B) = (R*inv(T))*Z' 44: * 45: * where inv(B) denotes the inverse of the matrix B, and Z' denotes the 46: * conjugate transpose of the matrix Z. 47: * 48: * Arguments 49: * ========= 50: * 51: * M (input) INTEGER 52: * The number of rows of the matrix A. M >= 0. 53: * 54: * P (input) INTEGER 55: * The number of rows of the matrix B. P >= 0. 56: * 57: * N (input) INTEGER 58: * The number of columns of the matrices A and B. N >= 0. 59: * 60: * A (input/output) COMPLEX*16 array, dimension (LDA,N) 61: * On entry, the M-by-N matrix A. 62: * On exit, if M <= N, the upper triangle of the subarray 63: * A(1:M,N-M+1:N) contains the M-by-M upper triangular matrix R; 64: * if M > N, the elements on and above the (M-N)-th subdiagonal 65: * contain the M-by-N upper trapezoidal matrix R; the remaining 66: * elements, with the array TAUA, represent the unitary 67: * matrix Q as a product of elementary reflectors (see Further 68: * Details). 69: * 70: * LDA (input) INTEGER 71: * The leading dimension of the array A. LDA >= max(1,M). 72: * 73: * TAUA (output) COMPLEX*16 array, dimension (min(M,N)) 74: * The scalar factors of the elementary reflectors which 75: * represent the unitary matrix Q (see Further Details). 76: * 77: * B (input/output) COMPLEX*16 array, dimension (LDB,N) 78: * On entry, the P-by-N matrix B. 79: * On exit, the elements on and above the diagonal of the array 80: * contain the min(P,N)-by-N upper trapezoidal matrix T (T is 81: * upper triangular if P >= N); the elements below the diagonal, 82: * with the array TAUB, represent the unitary matrix Z as a 83: * product of elementary reflectors (see Further Details). 84: * 85: * LDB (input) INTEGER 86: * The leading dimension of the array B. LDB >= max(1,P). 87: * 88: * TAUB (output) COMPLEX*16 array, dimension (min(P,N)) 89: * The scalar factors of the elementary reflectors which 90: * represent the unitary matrix Z (see Further Details). 91: * 92: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 93: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 94: * 95: * LWORK (input) INTEGER 96: * The dimension of the array WORK. LWORK >= max(1,N,M,P). 97: * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3), 98: * where NB1 is the optimal blocksize for the RQ factorization 99: * of an M-by-N matrix, NB2 is the optimal blocksize for the 100: * QR factorization of a P-by-N matrix, and NB3 is the optimal 101: * blocksize for a call of ZUNMRQ. 102: * 103: * If LWORK = -1, then a workspace query is assumed; the routine 104: * only calculates the optimal size of the WORK array, returns 105: * this value as the first entry of the WORK array, and no error 106: * message related to LWORK is issued by XERBLA. 107: * 108: * INFO (output) INTEGER 109: * = 0: successful exit 110: * < 0: if INFO=-i, the i-th argument had an illegal value. 111: * 112: * Further Details 113: * =============== 114: * 115: * The matrix Q is represented as a product of elementary reflectors 116: * 117: * Q = H(1) H(2) . . . H(k), where k = min(m,n). 118: * 119: * Each H(i) has the form 120: * 121: * H(i) = I - taua * v * v' 122: * 123: * where taua is a complex scalar, and v is a complex vector with 124: * v(n-k+i+1:n) = 0 and v(n-k+i) = 1; v(1:n-k+i-1) is stored on exit in 125: * A(m-k+i,1:n-k+i-1), and taua in TAUA(i). 126: * To form Q explicitly, use LAPACK subroutine ZUNGRQ. 127: * To use Q to update another matrix, use LAPACK subroutine ZUNMRQ. 128: * 129: * The matrix Z is represented as a product of elementary reflectors 130: * 131: * Z = H(1) H(2) . . . H(k), where k = min(p,n). 132: * 133: * Each H(i) has the form 134: * 135: * H(i) = I - taub * v * v' 136: * 137: * where taub is a complex scalar, and v is a complex vector with 138: * v(1:i-1) = 0 and v(i) = 1; v(i+1:p) is stored on exit in B(i+1:p,i), 139: * and taub in TAUB(i). 140: * To form Z explicitly, use LAPACK subroutine ZUNGQR. 141: * To use Z to update another matrix, use LAPACK subroutine ZUNMQR. 142: * 143: * ===================================================================== 144: * 145: * .. Local Scalars .. 146: LOGICAL LQUERY 147: INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3 148: * .. 149: * .. External Subroutines .. 150: EXTERNAL XERBLA, ZGEQRF, ZGERQF, ZUNMRQ 151: * .. 152: * .. External Functions .. 153: INTEGER ILAENV 154: EXTERNAL ILAENV 155: * .. 156: * .. Intrinsic Functions .. 157: INTRINSIC INT, MAX, MIN 158: * .. 159: * .. Executable Statements .. 160: * 161: * Test the input parameters 162: * 163: INFO = 0 164: NB1 = ILAENV( 1, 'ZGERQF', ' ', M, N, -1, -1 ) 165: NB2 = ILAENV( 1, 'ZGEQRF', ' ', P, N, -1, -1 ) 166: NB3 = ILAENV( 1, 'ZUNMRQ', ' ', M, N, P, -1 ) 167: NB = MAX( NB1, NB2, NB3 ) 168: LWKOPT = MAX( N, M, P )*NB 169: WORK( 1 ) = LWKOPT 170: LQUERY = ( LWORK.EQ.-1 ) 171: IF( M.LT.0 ) THEN 172: INFO = -1 173: ELSE IF( P.LT.0 ) THEN 174: INFO = -2 175: ELSE IF( N.LT.0 ) THEN 176: INFO = -3 177: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 178: INFO = -5 179: ELSE IF( LDB.LT.MAX( 1, P ) ) THEN 180: INFO = -8 181: ELSE IF( LWORK.LT.MAX( 1, M, P, N ) .AND. .NOT.LQUERY ) THEN 182: INFO = -11 183: END IF 184: IF( INFO.NE.0 ) THEN 185: CALL XERBLA( 'ZGGRQF', -INFO ) 186: RETURN 187: ELSE IF( LQUERY ) THEN 188: RETURN 189: END IF 190: * 191: * RQ factorization of M-by-N matrix A: A = R*Q 192: * 193: CALL ZGERQF( M, N, A, LDA, TAUA, WORK, LWORK, INFO ) 194: LOPT = WORK( 1 ) 195: * 196: * Update B := B*Q' 197: * 198: CALL ZUNMRQ( 'Right', 'Conjugate Transpose', P, N, MIN( M, N ), 199: $ A( MAX( 1, M-N+1 ), 1 ), LDA, TAUA, B, LDB, WORK, 200: $ LWORK, INFO ) 201: LOPT = MAX( LOPT, INT( WORK( 1 ) ) ) 202: * 203: * QR factorization of P-by-N matrix B: B = Z*T 204: * 205: CALL ZGEQRF( P, N, B, LDB, TAUB, WORK, LWORK, INFO ) 206: WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) ) 207: * 208: RETURN 209: * 210: * End of ZGGRQF 211: * 212: END