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