1 SUBROUTINE DORMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, 2 $ LDC, WORK, LWORK, INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER SIDE, TRANS, VECT 11 INTEGER INFO, K, LDA, LDC, LWORK, M, N 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) 15* .. 16* 17* Purpose 18* ======= 19* 20* If VECT = 'Q', DORMBR overwrites the general real M-by-N matrix C 21* with 22* SIDE = 'L' SIDE = 'R' 23* TRANS = 'N': Q * C C * Q 24* TRANS = 'T': Q**T * C C * Q**T 25* 26* If VECT = 'P', DORMBR overwrites the general real M-by-N matrix C 27* with 28* SIDE = 'L' SIDE = 'R' 29* TRANS = 'N': P * C C * P 30* TRANS = 'T': P**T * C C * P**T 31* 32* Here Q and P**T are the orthogonal matrices determined by DGEBRD when 33* reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and 34* P**T are defined as products of elementary reflectors H(i) and G(i) 35* respectively. 36* 37* Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the 38* order of the orthogonal matrix Q or P**T that is applied. 39* 40* If VECT = 'Q', A is assumed to have been an NQ-by-K matrix: 41* if nq >= k, Q = H(1) H(2) . . . H(k); 42* if nq < k, Q = H(1) H(2) . . . H(nq-1). 43* 44* If VECT = 'P', A is assumed to have been a K-by-NQ matrix: 45* if k < nq, P = G(1) G(2) . . . G(k); 46* if k >= nq, P = G(1) G(2) . . . G(nq-1). 47* 48* Arguments 49* ========= 50* 51* VECT (input) CHARACTER*1 52* = 'Q': apply Q or Q**T; 53* = 'P': apply P or P**T. 54* 55* SIDE (input) CHARACTER*1 56* = 'L': apply Q, Q**T, P or P**T from the Left; 57* = 'R': apply Q, Q**T, P or P**T from the Right. 58* 59* TRANS (input) CHARACTER*1 60* = 'N': No transpose, apply Q or P; 61* = 'T': Transpose, apply Q**T or P**T. 62* 63* M (input) INTEGER 64* The number of rows of the matrix C. M >= 0. 65* 66* N (input) INTEGER 67* The number of columns of the matrix C. N >= 0. 68* 69* K (input) INTEGER 70* If VECT = 'Q', the number of columns in the original 71* matrix reduced by DGEBRD. 72* If VECT = 'P', the number of rows in the original 73* matrix reduced by DGEBRD. 74* K >= 0. 75* 76* A (input) DOUBLE PRECISION array, dimension 77* (LDA,min(nq,K)) if VECT = 'Q' 78* (LDA,nq) if VECT = 'P' 79* The vectors which define the elementary reflectors H(i) and 80* G(i), whose products determine the matrices Q and P, as 81* returned by DGEBRD. 82* 83* LDA (input) INTEGER 84* The leading dimension of the array A. 85* If VECT = 'Q', LDA >= max(1,nq); 86* if VECT = 'P', LDA >= max(1,min(nq,K)). 87* 88* TAU (input) DOUBLE PRECISION array, dimension (min(nq,K)) 89* TAU(i) must contain the scalar factor of the elementary 90* reflector H(i) or G(i) which determines Q or P, as returned 91* by DGEBRD in the array argument TAUQ or TAUP. 92* 93* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 94* On entry, the M-by-N matrix C. 95* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q 96* or P*C or P**T*C or C*P or C*P**T. 97* 98* LDC (input) INTEGER 99* The leading dimension of the array C. LDC >= max(1,M). 100* 101* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 102* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 103* 104* LWORK (input) INTEGER 105* The dimension of the array WORK. 106* If SIDE = 'L', LWORK >= max(1,N); 107* if SIDE = 'R', LWORK >= max(1,M). 108* For optimum performance LWORK >= N*NB if SIDE = 'L', and 109* LWORK >= M*NB if SIDE = 'R', where NB is the optimal 110* blocksize. 111* 112* If LWORK = -1, then a workspace query is assumed; the routine 113* only calculates the optimal size of the WORK array, returns 114* this value as the first entry of the WORK array, and no error 115* message related to LWORK is issued by XERBLA. 116* 117* INFO (output) INTEGER 118* = 0: successful exit 119* < 0: if INFO = -i, the i-th argument had an illegal value 120* 121* ===================================================================== 122* 123* .. Local Scalars .. 124 LOGICAL APPLYQ, LEFT, LQUERY, NOTRAN 125 CHARACTER TRANST 126 INTEGER I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW 127* .. 128* .. External Functions .. 129 LOGICAL LSAME 130 INTEGER ILAENV 131 EXTERNAL LSAME, ILAENV 132* .. 133* .. External Subroutines .. 134 EXTERNAL DORMLQ, DORMQR, XERBLA 135* .. 136* .. Intrinsic Functions .. 137 INTRINSIC MAX, MIN 138* .. 139* .. Executable Statements .. 140* 141* Test the input arguments 142* 143 INFO = 0 144 APPLYQ = LSAME( VECT, 'Q' ) 145 LEFT = LSAME( SIDE, 'L' ) 146 NOTRAN = LSAME( TRANS, 'N' ) 147 LQUERY = ( LWORK.EQ.-1 ) 148* 149* NQ is the order of Q or P and NW is the minimum dimension of WORK 150* 151 IF( LEFT ) THEN 152 NQ = M 153 NW = N 154 ELSE 155 NQ = N 156 NW = M 157 END IF 158 IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN 159 INFO = -1 160 ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 161 INFO = -2 162 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN 163 INFO = -3 164 ELSE IF( M.LT.0 ) THEN 165 INFO = -4 166 ELSE IF( N.LT.0 ) THEN 167 INFO = -5 168 ELSE IF( K.LT.0 ) THEN 169 INFO = -6 170 ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR. 171 $ ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) ) 172 $ THEN 173 INFO = -8 174 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 175 INFO = -11 176 ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN 177 INFO = -13 178 END IF 179* 180 IF( INFO.EQ.0 ) THEN 181 IF( APPLYQ ) THEN 182 IF( LEFT ) THEN 183 NB = ILAENV( 1, 'DORMQR', SIDE // TRANS, M-1, N, M-1, 184 $ -1 ) 185 ELSE 186 NB = ILAENV( 1, 'DORMQR', SIDE // TRANS, M, N-1, N-1, 187 $ -1 ) 188 END IF 189 ELSE 190 IF( LEFT ) THEN 191 NB = ILAENV( 1, 'DORMLQ', SIDE // TRANS, M-1, N, M-1, 192 $ -1 ) 193 ELSE 194 NB = ILAENV( 1, 'DORMLQ', SIDE // TRANS, M, N-1, N-1, 195 $ -1 ) 196 END IF 197 END IF 198 LWKOPT = MAX( 1, NW )*NB 199 WORK( 1 ) = LWKOPT 200 END IF 201* 202 IF( INFO.NE.0 ) THEN 203 CALL XERBLA( 'DORMBR', -INFO ) 204 RETURN 205 ELSE IF( LQUERY ) THEN 206 RETURN 207 END IF 208* 209* Quick return if possible 210* 211 WORK( 1 ) = 1 212 IF( M.EQ.0 .OR. N.EQ.0 ) 213 $ RETURN 214* 215 IF( APPLYQ ) THEN 216* 217* Apply Q 218* 219 IF( NQ.GE.K ) THEN 220* 221* Q was determined by a call to DGEBRD with nq >= k 222* 223 CALL DORMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, 224 $ WORK, LWORK, IINFO ) 225 ELSE IF( NQ.GT.1 ) THEN 226* 227* Q was determined by a call to DGEBRD with nq < k 228* 229 IF( LEFT ) THEN 230 MI = M - 1 231 NI = N 232 I1 = 2 233 I2 = 1 234 ELSE 235 MI = M 236 NI = N - 1 237 I1 = 1 238 I2 = 2 239 END IF 240 CALL DORMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU, 241 $ C( I1, I2 ), LDC, WORK, LWORK, IINFO ) 242 END IF 243 ELSE 244* 245* Apply P 246* 247 IF( NOTRAN ) THEN 248 TRANST = 'T' 249 ELSE 250 TRANST = 'N' 251 END IF 252 IF( NQ.GT.K ) THEN 253* 254* P was determined by a call to DGEBRD with nq > k 255* 256 CALL DORMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC, 257 $ WORK, LWORK, IINFO ) 258 ELSE IF( NQ.GT.1 ) THEN 259* 260* P was determined by a call to DGEBRD with nq <= k 261* 262 IF( LEFT ) THEN 263 MI = M - 1 264 NI = N 265 I1 = 2 266 I2 = 1 267 ELSE 268 MI = M 269 NI = N - 1 270 I1 = 1 271 I2 = 2 272 END IF 273 CALL DORMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA, 274 $ TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO ) 275 END IF 276 END IF 277 WORK( 1 ) = LWKOPT 278 RETURN 279* 280* End of DORMBR 281* 282 END 283