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