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