00001 SUBROUTINE ZGGQRF( N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, 00002 $ LWORK, INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 INTEGER INFO, LDA, LDB, LWORK, M, N, P 00011 * .. 00012 * .. Array Arguments .. 00013 COMPLEX*16 A( LDA, * ), B( LDB, * ), TAUA( * ), TAUB( * ), 00014 $ WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * ZGGQRF computes a generalized QR factorization of an N-by-M matrix A 00021 * and an N-by-P matrix B: 00022 * 00023 * A = Q*R, B = Q*T*Z, 00024 * 00025 * where Q is an N-by-N unitary matrix, Z is a P-by-P unitary matrix, 00026 * and R and T assume one of the forms: 00027 * 00028 * if N >= M, R = ( R11 ) M , or if N < M, R = ( R11 R12 ) N, 00029 * ( 0 ) N-M N M-N 00030 * M 00031 * 00032 * where R11 is upper triangular, and 00033 * 00034 * if N <= P, T = ( 0 T12 ) N, or if N > P, T = ( T11 ) N-P, 00035 * P-N N ( T21 ) P 00036 * P 00037 * 00038 * where T12 or T21 is upper triangular. 00039 * 00040 * In particular, if B is square and nonsingular, the GQR factorization 00041 * of A and B implicitly gives the QR factorization of inv(B)*A: 00042 * 00043 * inv(B)*A = Z'*(inv(T)*R) 00044 * 00045 * where inv(B) denotes the inverse of the matrix B, and Z' denotes the 00046 * conjugate transpose of matrix Z. 00047 * 00048 * Arguments 00049 * ========= 00050 * 00051 * N (input) INTEGER 00052 * The number of rows of the matrices A and B. N >= 0. 00053 * 00054 * M (input) INTEGER 00055 * The number of columns of the matrix A. M >= 0. 00056 * 00057 * P (input) INTEGER 00058 * The number of columns of the matrix B. P >= 0. 00059 * 00060 * A (input/output) COMPLEX*16 array, dimension (LDA,M) 00061 * On entry, the N-by-M matrix A. 00062 * On exit, the elements on and above the diagonal of the array 00063 * contain the min(N,M)-by-M upper trapezoidal matrix R (R is 00064 * upper triangular if N >= M); the elements below the diagonal, 00065 * with the array TAUA, represent the unitary matrix Q as a 00066 * product of min(N,M) elementary reflectors (see Further 00067 * Details). 00068 * 00069 * LDA (input) INTEGER 00070 * The leading dimension of the array A. LDA >= max(1,N). 00071 * 00072 * TAUA (output) COMPLEX*16 array, dimension (min(N,M)) 00073 * The scalar factors of the elementary reflectors which 00074 * represent the unitary matrix Q (see Further Details). 00075 * 00076 * B (input/output) COMPLEX*16 array, dimension (LDB,P) 00077 * On entry, the N-by-P matrix B. 00078 * On exit, if N <= P, the upper triangle of the subarray 00079 * B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T; 00080 * if N > P, the elements on and above the (N-P)-th subdiagonal 00081 * contain the N-by-P upper trapezoidal matrix T; the remaining 00082 * elements, with the array TAUB, represent the unitary 00083 * matrix Z as a product of elementary reflectors (see Further 00084 * Details). 00085 * 00086 * LDB (input) INTEGER 00087 * The leading dimension of the array B. LDB >= max(1,N). 00088 * 00089 * TAUB (output) COMPLEX*16 array, dimension (min(N,P)) 00090 * The scalar factors of the elementary reflectors which 00091 * represent the unitary matrix Z (see Further Details). 00092 * 00093 * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK)) 00094 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00095 * 00096 * LWORK (input) INTEGER 00097 * The dimension of the array WORK. LWORK >= max(1,N,M,P). 00098 * For optimum performance LWORK >= max(N,M,P)*max(NB1,NB2,NB3), 00099 * where NB1 is the optimal blocksize for the QR factorization 00100 * of an N-by-M matrix, NB2 is the optimal blocksize for the 00101 * RQ factorization of an N-by-P matrix, and NB3 is the optimal 00102 * blocksize for a call of ZUNMQR. 00103 * 00104 * If LWORK = -1, then a workspace query is assumed; the routine 00105 * only calculates the optimal size of the WORK array, returns 00106 * this value as the first entry of the WORK array, and no error 00107 * message related to LWORK is issued by XERBLA. 00108 * 00109 * INFO (output) INTEGER 00110 * = 0: successful exit 00111 * < 0: if INFO = -i, the i-th argument had an illegal value. 00112 * 00113 * Further Details 00114 * =============== 00115 * 00116 * The matrix Q is represented as a product of elementary reflectors 00117 * 00118 * Q = H(1) H(2) . . . H(k), where k = min(n,m). 00119 * 00120 * Each H(i) has the form 00121 * 00122 * H(i) = I - taua * v * v' 00123 * 00124 * where taua is a complex scalar, and v is a complex vector with 00125 * v(1:i-1) = 0 and v(i) = 1; v(i+1:n) is stored on exit in A(i+1:n,i), 00126 * and taua in TAUA(i). 00127 * To form Q explicitly, use LAPACK subroutine ZUNGQR. 00128 * To use Q to update another matrix, use LAPACK subroutine ZUNMQR. 00129 * 00130 * The matrix Z is represented as a product of elementary reflectors 00131 * 00132 * Z = H(1) H(2) . . . H(k), where k = min(n,p). 00133 * 00134 * Each H(i) has the form 00135 * 00136 * H(i) = I - taub * v * v' 00137 * 00138 * where taub is a complex scalar, and v is a complex vector with 00139 * v(p-k+i+1:p) = 0 and v(p-k+i) = 1; v(1:p-k+i-1) is stored on exit in 00140 * B(n-k+i,1:p-k+i-1), and taub in TAUB(i). 00141 * To form Z explicitly, use LAPACK subroutine ZUNGRQ. 00142 * To use Z to update another matrix, use LAPACK subroutine ZUNMRQ. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. Local Scalars .. 00147 LOGICAL LQUERY 00148 INTEGER LOPT, LWKOPT, NB, NB1, NB2, NB3 00149 * .. 00150 * .. External Subroutines .. 00151 EXTERNAL XERBLA, ZGEQRF, ZGERQF, ZUNMQR 00152 * .. 00153 * .. External Functions .. 00154 INTEGER ILAENV 00155 EXTERNAL ILAENV 00156 * .. 00157 * .. Intrinsic Functions .. 00158 INTRINSIC INT, MAX, MIN 00159 * .. 00160 * .. Executable Statements .. 00161 * 00162 * Test the input parameters 00163 * 00164 INFO = 0 00165 NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, M, -1, -1 ) 00166 NB2 = ILAENV( 1, 'ZGERQF', ' ', N, P, -1, -1 ) 00167 NB3 = ILAENV( 1, 'ZUNMQR', ' ', N, M, P, -1 ) 00168 NB = MAX( NB1, NB2, NB3 ) 00169 LWKOPT = MAX( N, M, P )*NB 00170 WORK( 1 ) = LWKOPT 00171 LQUERY = ( LWORK.EQ.-1 ) 00172 IF( N.LT.0 ) THEN 00173 INFO = -1 00174 ELSE IF( M.LT.0 ) THEN 00175 INFO = -2 00176 ELSE IF( P.LT.0 ) THEN 00177 INFO = -3 00178 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00179 INFO = -5 00180 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00181 INFO = -8 00182 ELSE IF( LWORK.LT.MAX( 1, N, M, P ) .AND. .NOT.LQUERY ) THEN 00183 INFO = -11 00184 END IF 00185 IF( INFO.NE.0 ) THEN 00186 CALL XERBLA( 'ZGGQRF', -INFO ) 00187 RETURN 00188 ELSE IF( LQUERY ) THEN 00189 RETURN 00190 END IF 00191 * 00192 * QR factorization of N-by-M matrix A: A = Q*R 00193 * 00194 CALL ZGEQRF( N, M, A, LDA, TAUA, WORK, LWORK, INFO ) 00195 LOPT = WORK( 1 ) 00196 * 00197 * Update B := Q'*B. 00198 * 00199 CALL ZUNMQR( 'Left', 'Conjugate Transpose', N, P, MIN( N, M ), A, 00200 $ LDA, TAUA, B, LDB, WORK, LWORK, INFO ) 00201 LOPT = MAX( LOPT, INT( WORK( 1 ) ) ) 00202 * 00203 * RQ factorization of N-by-P matrix B: B = T*Z. 00204 * 00205 CALL ZGERQF( N, P, B, LDB, TAUB, WORK, LWORK, INFO ) 00206 WORK( 1 ) = MAX( LOPT, INT( WORK( 1 ) ) ) 00207 * 00208 RETURN 00209 * 00210 * End of ZGGQRF 00211 * 00212 END
1.7.2