 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, 00002 $ LDV, T, LDT, C, LDC, WORK, LDWORK ) 00003 * 00004 * -- LAPACK routine (version 3.3.1) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * -- April 2011 -- 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER DIRECT, SIDE, STOREV, TRANS 00011 INTEGER K, L, LDC, LDT, LDV, LDWORK, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 REAL C( LDC, * ), T( LDT, * ), V( LDV, * ), 00015 $ WORK( LDWORK, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * SLARZB applies a real block reflector H or its transpose H**T to 00022 * a real distributed M-by-N C from the left or the right. 00023 * 00024 * Currently, only STOREV = 'R' and DIRECT = 'B' are supported. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * SIDE (input) CHARACTER*1 00030 * = 'L': apply H or H**T from the Left 00031 * = 'R': apply H or H**T from the Right 00032 * 00033 * TRANS (input) CHARACTER*1 00034 * = 'N': apply H (No transpose) 00035 * = 'C': apply H**T (Transpose) 00036 * 00037 * DIRECT (input) CHARACTER*1 00038 * Indicates how H is formed from a product of elementary 00039 * reflectors 00040 * = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet) 00041 * = 'B': H = H(k) . . . H(2) H(1) (Backward) 00042 * 00043 * STOREV (input) CHARACTER*1 00044 * Indicates how the vectors which define the elementary 00045 * reflectors are stored: 00046 * = 'C': Columnwise (not supported yet) 00047 * = 'R': Rowwise 00048 * 00049 * M (input) INTEGER 00050 * The number of rows of the matrix C. 00051 * 00052 * N (input) INTEGER 00053 * The number of columns of the matrix C. 00054 * 00055 * K (input) INTEGER 00056 * The order of the matrix T (= the number of elementary 00057 * reflectors whose product defines the block reflector). 00058 * 00059 * L (input) INTEGER 00060 * The number of columns of the matrix V containing the 00061 * meaningful part of the Householder reflectors. 00062 * If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0. 00063 * 00064 * V (input) REAL array, dimension (LDV,NV). 00065 * If STOREV = 'C', NV = K; if STOREV = 'R', NV = L. 00066 * 00067 * LDV (input) INTEGER 00068 * The leading dimension of the array V. 00069 * If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K. 00070 * 00071 * T (input) REAL array, dimension (LDT,K) 00072 * The triangular K-by-K matrix T in the representation of the 00073 * block reflector. 00074 * 00075 * LDT (input) INTEGER 00076 * The leading dimension of the array T. LDT >= K. 00077 * 00078 * C (input/output) REAL array, dimension (LDC,N) 00079 * On entry, the M-by-N matrix C. 00080 * On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T. 00081 * 00082 * LDC (input) INTEGER 00083 * The leading dimension of the array C. LDC >= max(1,M). 00084 * 00085 * WORK (workspace) REAL array, dimension (LDWORK,K) 00086 * 00087 * LDWORK (input) INTEGER 00088 * The leading dimension of the array WORK. 00089 * If SIDE = 'L', LDWORK >= max(1,N); 00090 * if SIDE = 'R', LDWORK >= max(1,M). 00091 * 00092 * Further Details 00093 * =============== 00094 * 00095 * Based on contributions by 00096 * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA 00097 * 00098 * ===================================================================== 00099 * 00100 * .. Parameters .. 00101 REAL ONE 00102 PARAMETER ( ONE = 1.0E+0 ) 00103 * .. 00104 * .. Local Scalars .. 00105 CHARACTER TRANST 00106 INTEGER I, INFO, J 00107 * .. 00108 * .. External Functions .. 00109 LOGICAL LSAME 00110 EXTERNAL LSAME 00111 * .. 00112 * .. External Subroutines .. 00113 EXTERNAL SCOPY, SGEMM, STRMM, XERBLA 00114 * .. 00115 * .. Executable Statements .. 00116 * 00117 * Quick return if possible 00118 * 00119 IF( M.LE.0 .OR. N.LE.0 ) 00120 $ RETURN 00121 * 00122 * Check for currently supported options 00123 * 00124 INFO = 0 00125 IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN 00126 INFO = -3 00127 ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN 00128 INFO = -4 00129 END IF 00130 IF( INFO.NE.0 ) THEN 00131 CALL XERBLA( 'SLARZB', -INFO ) 00132 RETURN 00133 END IF 00134 * 00135 IF( LSAME( TRANS, 'N' ) ) THEN 00136 TRANST = 'T' 00137 ELSE 00138 TRANST = 'N' 00139 END IF 00140 * 00141 IF( LSAME( SIDE, 'L' ) ) THEN 00142 * 00143 * Form H * C or H**T * C 00144 * 00145 * W( 1:n, 1:k ) = C( 1:k, 1:n )**T 00146 * 00147 DO 10 J = 1, K 00148 CALL SCOPY( N, C( J, 1 ), LDC, WORK( 1, J ), 1 ) 00149 10 CONTINUE 00150 * 00151 * W( 1:n, 1:k ) = W( 1:n, 1:k ) + ... 00152 * C( m-l+1:m, 1:n )**T * V( 1:k, 1:l )**T 00153 * 00154 IF( L.GT.0 ) 00155 $ CALL SGEMM( 'Transpose', 'Transpose', N, K, L, ONE, 00156 $ C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, LDWORK ) 00157 * 00158 * W( 1:n, 1:k ) = W( 1:n, 1:k ) * T**T or W( 1:m, 1:k ) * T 00159 * 00160 CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T, 00161 $ LDT, WORK, LDWORK ) 00162 * 00163 * C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )**T 00164 * 00165 DO 30 J = 1, N 00166 DO 20 I = 1, K 00167 C( I, J ) = C( I, J ) - WORK( J, I ) 00168 20 CONTINUE 00169 30 CONTINUE 00170 * 00171 * C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ... 00172 * V( 1:k, 1:l )**T * W( 1:n, 1:k )**T 00173 * 00174 IF( L.GT.0 ) 00175 $ CALL SGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV, 00176 $ WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC ) 00177 * 00178 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00179 * 00180 * Form C * H or C * H**T 00181 * 00182 * W( 1:m, 1:k ) = C( 1:m, 1:k ) 00183 * 00184 DO 40 J = 1, K 00185 CALL SCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 ) 00186 40 CONTINUE 00187 * 00188 * W( 1:m, 1:k ) = W( 1:m, 1:k ) + ... 00189 * C( 1:m, n-l+1:n ) * V( 1:k, 1:l )**T 00190 * 00191 IF( L.GT.0 ) 00192 $ CALL SGEMM( 'No transpose', 'Transpose', M, K, L, ONE, 00193 $ C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK ) 00194 * 00195 * W( 1:m, 1:k ) = W( 1:m, 1:k ) * T or W( 1:m, 1:k ) * T**T 00196 * 00197 CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T, 00198 $ LDT, WORK, LDWORK ) 00199 * 00200 * C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k ) 00201 * 00202 DO 60 J = 1, K 00203 DO 50 I = 1, M 00204 C( I, J ) = C( I, J ) - WORK( I, J ) 00205 50 CONTINUE 00206 60 CONTINUE 00207 * 00208 * C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ... 00209 * W( 1:m, 1:k ) * V( 1:k, 1:l ) 00210 * 00211 IF( L.GT.0 ) 00212 $ CALL SGEMM( 'No transpose', 'No transpose', M, L, K, -ONE, 00213 $ WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC ) 00214 * 00215 END IF 00216 * 00217 RETURN 00218 * 00219 * End of SLARZB 00220 * 00221 END
1.7.3 
