 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, 00002 $ INFO ) 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 SIDE, TRANS, UPLO 00011 INTEGER INFO, LDC, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 COMPLEX AP( * ), C( LDC, * ), TAU( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * CUPMTR overwrites the general complex M-by-N matrix C with 00021 * 00022 * SIDE = 'L' SIDE = 'R' 00023 * TRANS = 'N': Q * C C * Q 00024 * TRANS = 'C': Q**H * C C * Q**H 00025 * 00026 * where Q is a complex unitary matrix of order nq, with nq = m if 00027 * SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of 00028 * nq-1 elementary reflectors, as returned by CHPTRD using packed 00029 * storage: 00030 * 00031 * if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1); 00032 * 00033 * if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1). 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * SIDE (input) CHARACTER*1 00039 * = 'L': apply Q or Q**H from the Left; 00040 * = 'R': apply Q or Q**H from the Right. 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': Upper triangular packed storage used in previous 00044 * call to CHPTRD; 00045 * = 'L': Lower triangular packed storage used in previous 00046 * call to CHPTRD. 00047 * 00048 * TRANS (input) CHARACTER*1 00049 * = 'N': No transpose, apply Q; 00050 * = 'C': Conjugate transpose, apply Q**H. 00051 * 00052 * M (input) INTEGER 00053 * The number of rows of the matrix C. M >= 0. 00054 * 00055 * N (input) INTEGER 00056 * The number of columns of the matrix C. N >= 0. 00057 * 00058 * AP (input) COMPLEX array, dimension 00059 * (M*(M+1)/2) if SIDE = 'L' 00060 * (N*(N+1)/2) if SIDE = 'R' 00061 * The vectors which define the elementary reflectors, as 00062 * returned by CHPTRD. AP is modified by the routine but 00063 * restored on exit. 00064 * 00065 * TAU (input) COMPLEX array, dimension (M-1) if SIDE = 'L' 00066 * or (N-1) if SIDE = 'R' 00067 * TAU(i) must contain the scalar factor of the elementary 00068 * reflector H(i), as returned by CHPTRD. 00069 * 00070 * C (input/output) COMPLEX array, dimension (LDC,N) 00071 * On entry, the M-by-N matrix C. 00072 * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q. 00073 * 00074 * LDC (input) INTEGER 00075 * The leading dimension of the array C. LDC >= max(1,M). 00076 * 00077 * WORK (workspace) COMPLEX array, dimension 00078 * (N) if SIDE = 'L' 00079 * (M) if SIDE = 'R' 00080 * 00081 * INFO (output) INTEGER 00082 * = 0: successful exit 00083 * < 0: if INFO = -i, the i-th argument had an illegal value 00084 * 00085 * ===================================================================== 00086 * 00087 * .. Parameters .. 00088 COMPLEX ONE 00089 PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) 00090 * .. 00091 * .. Local Scalars .. 00092 LOGICAL FORWRD, LEFT, NOTRAN, UPPER 00093 INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ 00094 COMPLEX AII, TAUI 00095 * .. 00096 * .. External Functions .. 00097 LOGICAL LSAME 00098 EXTERNAL LSAME 00099 * .. 00100 * .. External Subroutines .. 00101 EXTERNAL CLARF, XERBLA 00102 * .. 00103 * .. Intrinsic Functions .. 00104 INTRINSIC CONJG, MAX 00105 * .. 00106 * .. Executable Statements .. 00107 * 00108 * Test the input arguments 00109 * 00110 INFO = 0 00111 LEFT = LSAME( SIDE, 'L' ) 00112 NOTRAN = LSAME( TRANS, 'N' ) 00113 UPPER = LSAME( UPLO, 'U' ) 00114 * 00115 * NQ is the order of Q 00116 * 00117 IF( LEFT ) THEN 00118 NQ = M 00119 ELSE 00120 NQ = N 00121 END IF 00122 IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN 00123 INFO = -1 00124 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 00125 INFO = -2 00126 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN 00127 INFO = -3 00128 ELSE IF( M.LT.0 ) THEN 00129 INFO = -4 00130 ELSE IF( N.LT.0 ) THEN 00131 INFO = -5 00132 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 00133 INFO = -9 00134 END IF 00135 IF( INFO.NE.0 ) THEN 00136 CALL XERBLA( 'CUPMTR', -INFO ) 00137 RETURN 00138 END IF 00139 * 00140 * Quick return if possible 00141 * 00142 IF( M.EQ.0 .OR. N.EQ.0 ) 00143 $ RETURN 00144 * 00145 IF( UPPER ) THEN 00146 * 00147 * Q was determined by a call to CHPTRD with UPLO = 'U' 00148 * 00149 FORWRD = ( LEFT .AND. NOTRAN ) .OR. 00150 $ ( .NOT.LEFT .AND. .NOT.NOTRAN ) 00151 * 00152 IF( FORWRD ) THEN 00153 I1 = 1 00154 I2 = NQ - 1 00155 I3 = 1 00156 II = 2 00157 ELSE 00158 I1 = NQ - 1 00159 I2 = 1 00160 I3 = -1 00161 II = NQ*( NQ+1 ) / 2 - 1 00162 END IF 00163 * 00164 IF( LEFT ) THEN 00165 NI = N 00166 ELSE 00167 MI = M 00168 END IF 00169 * 00170 DO 10 I = I1, I2, I3 00171 IF( LEFT ) THEN 00172 * 00173 * H(i) or H(i)**H is applied to C(1:i,1:n) 00174 * 00175 MI = I 00176 ELSE 00177 * 00178 * H(i) or H(i)**H is applied to C(1:m,1:i) 00179 * 00180 NI = I 00181 END IF 00182 * 00183 * Apply H(i) or H(i)**H 00184 * 00185 IF( NOTRAN ) THEN 00186 TAUI = TAU( I ) 00187 ELSE 00188 TAUI = CONJG( TAU( I ) ) 00189 END IF 00190 AII = AP( II ) 00191 AP( II ) = ONE 00192 CALL CLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, LDC, 00193 $ WORK ) 00194 AP( II ) = AII 00195 * 00196 IF( FORWRD ) THEN 00197 II = II + I + 2 00198 ELSE 00199 II = II - I - 1 00200 END IF 00201 10 CONTINUE 00202 ELSE 00203 * 00204 * Q was determined by a call to CHPTRD with UPLO = 'L'. 00205 * 00206 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR. 00207 $ ( .NOT.LEFT .AND. NOTRAN ) 00208 * 00209 IF( FORWRD ) THEN 00210 I1 = 1 00211 I2 = NQ - 1 00212 I3 = 1 00213 II = 2 00214 ELSE 00215 I1 = NQ - 1 00216 I2 = 1 00217 I3 = -1 00218 II = NQ*( NQ+1 ) / 2 - 1 00219 END IF 00220 * 00221 IF( LEFT ) THEN 00222 NI = N 00223 JC = 1 00224 ELSE 00225 MI = M 00226 IC = 1 00227 END IF 00228 * 00229 DO 20 I = I1, I2, I3 00230 AII = AP( II ) 00231 AP( II ) = ONE 00232 IF( LEFT ) THEN 00233 * 00234 * H(i) or H(i)**H is applied to C(i+1:m,1:n) 00235 * 00236 MI = M - I 00237 IC = I + 1 00238 ELSE 00239 * 00240 * H(i) or H(i)**H is applied to C(1:m,i+1:n) 00241 * 00242 NI = N - I 00243 JC = I + 1 00244 END IF 00245 * 00246 * Apply H(i) or H(i)**H 00247 * 00248 IF( NOTRAN ) THEN 00249 TAUI = TAU( I ) 00250 ELSE 00251 TAUI = CONJG( TAU( I ) ) 00252 END IF 00253 CALL CLARF( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, JC ), 00254 $ LDC, WORK ) 00255 AP( II ) = AII 00256 * 00257 IF( FORWRD ) THEN 00258 II = II + NQ - I + 1 00259 ELSE 00260 II = II - NQ + I - 2 00261 END IF 00262 20 CONTINUE 00263 END IF 00264 RETURN 00265 * 00266 * End of CUPMTR 00267 * 00268 END
1.7.3 
