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