00001 SUBROUTINE CLARZB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, V, 00002 $ LDV, T, LDT, C, LDC, WORK, LDWORK ) 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 DIRECT, SIDE, STOREV, TRANS 00011 INTEGER K, L, LDC, LDT, LDV, LDWORK, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), 00015 $ WORK( LDWORK, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * CLARZB applies a complex block reflector H or its transpose H**H 00022 * to a complex 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' from the Left 00031 * = 'R': apply H or H' from the Right 00032 * 00033 * TRANS (input) CHARACTER*1 00034 * = 'N': apply H (No transpose) 00035 * = 'C': apply H' (Conjugate 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) COMPLEX 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) COMPLEX 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) COMPLEX array, dimension (LDC,N) 00079 * On entry, the M-by-N matrix C. 00080 * On exit, C is overwritten by H*C or H'*C or C*H or C*H'. 00081 * 00082 * LDC (input) INTEGER 00083 * The leading dimension of the array C. LDC >= max(1,M). 00084 * 00085 * WORK (workspace) COMPLEX 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 COMPLEX ONE 00102 PARAMETER ( ONE = ( 1.0E+0, 0.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 CCOPY, CGEMM, CLACGV, CTRMM, 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( 'CLARZB', -INFO ) 00132 RETURN 00133 END IF 00134 * 00135 IF( LSAME( TRANS, 'N' ) ) THEN 00136 TRANST = 'C' 00137 ELSE 00138 TRANST = 'N' 00139 END IF 00140 * 00141 IF( LSAME( SIDE, 'L' ) ) THEN 00142 * 00143 * Form H * C or H' * C 00144 * 00145 * W( 1:n, 1:k ) = conjg( C( 1:k, 1:n )' ) 00146 * 00147 DO 10 J = 1, K 00148 CALL CCOPY( 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 * conjg( C( m-l+1:m, 1:n )' ) * V( 1:k, 1:l )' 00153 * 00154 IF( L.GT.0 ) 00155 $ CALL CGEMM( 'Transpose', 'Conjugate transpose', N, K, L, 00156 $ ONE, C( M-L+1, 1 ), LDC, V, LDV, ONE, WORK, 00157 $ LDWORK ) 00158 * 00159 * W( 1:n, 1:k ) = W( 1:n, 1:k ) * T' or W( 1:m, 1:k ) * T 00160 * 00161 CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', N, K, ONE, T, 00162 $ LDT, WORK, LDWORK ) 00163 * 00164 * C( 1:k, 1:n ) = C( 1:k, 1:n ) - conjg( W( 1:n, 1:k )' ) 00165 * 00166 DO 30 J = 1, N 00167 DO 20 I = 1, K 00168 C( I, J ) = C( I, J ) - WORK( J, I ) 00169 20 CONTINUE 00170 30 CONTINUE 00171 * 00172 * C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ... 00173 * conjg( V( 1:k, 1:l )' ) * conjg( W( 1:n, 1:k )' ) 00174 * 00175 IF( L.GT.0 ) 00176 $ CALL CGEMM( 'Transpose', 'Transpose', L, N, K, -ONE, V, LDV, 00177 $ WORK, LDWORK, ONE, C( M-L+1, 1 ), LDC ) 00178 * 00179 ELSE IF( LSAME( SIDE, 'R' ) ) THEN 00180 * 00181 * Form C * H or C * H' 00182 * 00183 * W( 1:m, 1:k ) = C( 1:m, 1:k ) 00184 * 00185 DO 40 J = 1, K 00186 CALL CCOPY( M, C( 1, J ), 1, WORK( 1, J ), 1 ) 00187 40 CONTINUE 00188 * 00189 * W( 1:m, 1:k ) = W( 1:m, 1:k ) + ... 00190 * C( 1:m, n-l+1:n ) * conjg( V( 1:k, 1:l )' ) 00191 * 00192 IF( L.GT.0 ) 00193 $ CALL CGEMM( 'No transpose', 'Transpose', M, K, L, ONE, 00194 $ C( 1, N-L+1 ), LDC, V, LDV, ONE, WORK, LDWORK ) 00195 * 00196 * W( 1:m, 1:k ) = W( 1:m, 1:k ) * conjg( T ) or 00197 * W( 1:m, 1:k ) * conjg( T' ) 00198 * 00199 DO 50 J = 1, K 00200 CALL CLACGV( K-J+1, T( J, J ), 1 ) 00201 50 CONTINUE 00202 CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', M, K, ONE, T, 00203 $ LDT, WORK, LDWORK ) 00204 DO 60 J = 1, K 00205 CALL CLACGV( K-J+1, T( J, J ), 1 ) 00206 60 CONTINUE 00207 * 00208 * C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k ) 00209 * 00210 DO 80 J = 1, K 00211 DO 70 I = 1, M 00212 C( I, J ) = C( I, J ) - WORK( I, J ) 00213 70 CONTINUE 00214 80 CONTINUE 00215 * 00216 * C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ... 00217 * W( 1:m, 1:k ) * conjg( V( 1:k, 1:l ) ) 00218 * 00219 DO 90 J = 1, L 00220 CALL CLACGV( K, V( 1, J ), 1 ) 00221 90 CONTINUE 00222 IF( L.GT.0 ) 00223 $ CALL CGEMM( 'No transpose', 'No transpose', M, L, K, -ONE, 00224 $ WORK, LDWORK, V, LDV, ONE, C( 1, N-L+1 ), LDC ) 00225 DO 100 J = 1, L 00226 CALL CLACGV( K, V( 1, J ), 1 ) 00227 100 CONTINUE 00228 * 00229 END IF 00230 * 00231 RETURN 00232 * 00233 * End of CLARZB 00234 * 00235 END
1.7.2