LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CUNMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, 00002 $ LDC, WORK, 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 CHARACTER SIDE, TRANS, VECT 00011 INTEGER INFO, K, LDA, LDC, LWORK, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX A( LDA, * ), C( LDC, * ), TAU( * ), 00015 $ WORK( * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * If VECT = 'Q', CUNMBR overwrites the general complex M-by-N matrix C 00022 * with 00023 * SIDE = 'L' SIDE = 'R' 00024 * TRANS = 'N': Q * C C * Q 00025 * TRANS = 'C': Q**H * C C * Q**H 00026 * 00027 * If VECT = 'P', CUNMBR overwrites the general complex M-by-N matrix C 00028 * with 00029 * SIDE = 'L' SIDE = 'R' 00030 * TRANS = 'N': P * C C * P 00031 * TRANS = 'C': P**H * C C * P**H 00032 * 00033 * Here Q and P**H are the unitary matrices determined by CGEBRD when 00034 * reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q 00035 * and P**H are defined as products of elementary reflectors H(i) and 00036 * G(i) respectively. 00037 * 00038 * Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the 00039 * order of the unitary matrix Q or P**H that is applied. 00040 * 00041 * If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: 00042 * if nq >= k, Q = H(1) H(2) . . . H(k); 00043 * if nq < k, Q = H(1) H(2) . . . H(nq-1). 00044 * 00045 * If VECT = 'P', A is assumed to have been a K-by-NQ matrix: 00046 * if k < nq, P = G(1) G(2) . . . G(k); 00047 * if k >= nq, P = G(1) G(2) . . . G(nq-1). 00048 * 00049 * Arguments 00050 * ========= 00051 * 00052 * VECT (input) CHARACTER*1 00053 * = 'Q': apply Q or Q**H; 00054 * = 'P': apply P or P**H. 00055 * 00056 * SIDE (input) CHARACTER*1 00057 * = 'L': apply Q, Q**H, P or P**H from the Left; 00058 * = 'R': apply Q, Q**H, P or P**H from the Right. 00059 * 00060 * TRANS (input) CHARACTER*1 00061 * = 'N': No transpose, apply Q or P; 00062 * = 'C': Conjugate transpose, apply Q**H or P**H. 00063 * 00064 * M (input) INTEGER 00065 * The number of rows of the matrix C. M >= 0. 00066 * 00067 * N (input) INTEGER 00068 * The number of columns of the matrix C. N >= 0. 00069 * 00070 * K (input) INTEGER 00071 * If VECT = 'Q', the number of columns in the original 00072 * matrix reduced by CGEBRD. 00073 * If VECT = 'P', the number of rows in the original 00074 * matrix reduced by CGEBRD. 00075 * K >= 0. 00076 * 00077 * A (input) COMPLEX array, dimension 00078 * (LDA,min(nq,K)) if VECT = 'Q' 00079 * (LDA,nq) if VECT = 'P' 00080 * The vectors which define the elementary reflectors H(i) and 00081 * G(i), whose products determine the matrices Q and P, as 00082 * returned by CGEBRD. 00083 * 00084 * LDA (input) INTEGER 00085 * The leading dimension of the array A. 00086 * If VECT = 'Q', LDA >= max(1,nq); 00087 * if VECT = 'P', LDA >= max(1,min(nq,K)). 00088 * 00089 * TAU (input) COMPLEX array, dimension (min(nq,K)) 00090 * TAU(i) must contain the scalar factor of the elementary 00091 * reflector H(i) or G(i) which determines Q or P, as returned 00092 * by CGEBRD in the array argument TAUQ or TAUP. 00093 * 00094 * C (input/output) COMPLEX array, dimension (LDC,N) 00095 * On entry, the M-by-N matrix C. 00096 * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q 00097 * or P*C or P**H*C or C*P or C*P**H. 00098 * 00099 * LDC (input) INTEGER 00100 * The leading dimension of the array C. LDC >= max(1,M). 00101 * 00102 * WORK (workspace/output) COMPLEX array, dimension (MAX(1,LWORK)) 00103 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00104 * 00105 * LWORK (input) INTEGER 00106 * The dimension of the array WORK. 00107 * If SIDE = 'L', LWORK >= max(1,N); 00108 * if SIDE = 'R', LWORK >= max(1,M); 00109 * if N = 0 or M = 0, LWORK >= 1. 00110 * For optimum performance LWORK >= max(1,N*NB) if SIDE = 'L', 00111 * and LWORK >= max(1,M*NB) if SIDE = 'R', where NB is the 00112 * optimal blocksize. (NB = 0 if M = 0 or N = 0.) 00113 * 00114 * If LWORK = -1, then a workspace query is assumed; the routine 00115 * only calculates the optimal size of the WORK array, returns 00116 * this value as the first entry of the WORK array, and no error 00117 * message related to LWORK is issued by XERBLA. 00118 * 00119 * INFO (output) INTEGER 00120 * = 0: successful exit 00121 * < 0: if INFO = -i, the i-th argument had an illegal value 00122 * 00123 * ===================================================================== 00124 * 00125 * .. Local Scalars .. 00126 LOGICAL APPLYQ, LEFT, LQUERY, NOTRAN 00127 CHARACTER TRANST 00128 INTEGER I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW 00129 * .. 00130 * .. External Functions .. 00131 LOGICAL LSAME 00132 INTEGER ILAENV 00133 EXTERNAL ILAENV, LSAME 00134 * .. 00135 * .. External Subroutines .. 00136 EXTERNAL CUNMLQ, CUNMQR, XERBLA 00137 * .. 00138 * .. Intrinsic Functions .. 00139 INTRINSIC MAX, MIN 00140 * .. 00141 * .. Executable Statements .. 00142 * 00143 * Test the input arguments 00144 * 00145 INFO = 0 00146 APPLYQ = LSAME( VECT, 'Q' ) 00147 LEFT = LSAME( SIDE, 'L' ) 00148 NOTRAN = LSAME( TRANS, 'N' ) 00149 LQUERY = ( LWORK.EQ.-1 ) 00150 * 00151 * NQ is the order of Q or P and NW is the minimum dimension of WORK 00152 * 00153 IF( LEFT ) THEN 00154 NQ = M 00155 NW = N 00156 ELSE 00157 NQ = N 00158 NW = M 00159 END IF 00160 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00161 NW = 0 00162 END IF 00163 IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN 00164 INFO = -1 00165 ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 00166 INFO = -2 00167 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00168 INFO = -3 00169 ELSE IF( M.LT.0 ) THEN 00170 INFO = -4 00171 ELSE IF( N.LT.0 ) THEN 00172 INFO = -5 00173 ELSE IF( K.LT.0 ) THEN 00174 INFO = -6 00175 ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR. 00176 $ ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) ) 00177 $ THEN 00178 INFO = -8 00179 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00180 INFO = -11 00181 ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN 00182 INFO = -13 00183 END IF 00184 * 00185 IF( INFO.EQ.0 ) THEN 00186 IF( NW.GT.0 ) THEN 00187 IF( APPLYQ ) THEN 00188 IF( LEFT ) THEN 00189 NB = ILAENV( 1, 'CUNMQR', SIDE // TRANS, M-1, N, M-1, 00190 $ -1 ) 00191 ELSE 00192 NB = ILAENV( 1, 'CUNMQR', SIDE // TRANS, M, N-1, N-1, 00193 $ -1 ) 00194 END IF 00195 ELSE 00196 IF( LEFT ) THEN 00197 NB = ILAENV( 1, 'CUNMLQ', SIDE // TRANS, M-1, N, M-1, 00198 $ -1 ) 00199 ELSE 00200 NB = ILAENV( 1, 'CUNMLQ', SIDE // TRANS, M, N-1, N-1, 00201 $ -1 ) 00202 END IF 00203 END IF 00204 LWKOPT = MAX( 1, NW*NB ) 00205 ELSE 00206 LWKOPT = 1 00207 END IF 00208 WORK( 1 ) = LWKOPT 00209 END IF 00210 * 00211 IF( INFO.NE.0 ) THEN 00212 CALL XERBLA( 'CUNMBR', -INFO ) 00213 RETURN 00214 ELSE IF( LQUERY ) THEN 00215 RETURN 00216 END IF 00217 * 00218 * Quick return if possible 00219 * 00220 IF( M.EQ.0 .OR. N.EQ.0 ) 00221 $ RETURN 00222 * 00223 IF( APPLYQ ) THEN 00224 * 00225 * Apply Q 00226 * 00227 IF( NQ.GE.K ) THEN 00228 * 00229 * Q was determined by a call to CGEBRD with nq >= k 00230 * 00231 CALL CUNMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, 00232 $ WORK, LWORK, IINFO ) 00233 ELSE IF( NQ.GT.1 ) THEN 00234 * 00235 * Q was determined by a call to CGEBRD with nq < k 00236 * 00237 IF( LEFT ) THEN 00238 MI = M - 1 00239 NI = N 00240 I1 = 2 00241 I2 = 1 00242 ELSE 00243 MI = M 00244 NI = N - 1 00245 I1 = 1 00246 I2 = 2 00247 END IF 00248 CALL CUNMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU, 00249 $ C( I1, I2 ), LDC, WORK, LWORK, IINFO ) 00250 END IF 00251 ELSE 00252 * 00253 * Apply P 00254 * 00255 IF( NOTRAN ) THEN 00256 TRANST = 'C' 00257 ELSE 00258 TRANST = 'N' 00259 END IF 00260 IF( NQ.GT.K ) THEN 00261 * 00262 * P was determined by a call to CGEBRD with nq > k 00263 * 00264 CALL CUNMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC, 00265 $ WORK, LWORK, IINFO ) 00266 ELSE IF( NQ.GT.1 ) THEN 00267 * 00268 * P was determined by a call to CGEBRD with nq <= k 00269 * 00270 IF( LEFT ) THEN 00271 MI = M - 1 00272 NI = N 00273 I1 = 2 00274 I2 = 1 00275 ELSE 00276 MI = M 00277 NI = N - 1 00278 I1 = 1 00279 I2 = 2 00280 END IF 00281 CALL CUNMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA, 00282 $ TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO ) 00283 END IF 00284 END IF 00285 WORK( 1 ) = LWKOPT 00286 RETURN 00287 * 00288 * End of CUNMBR 00289 * 00290 END