|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DOPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, 00002 $ 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, UPLO 00011 INTEGER INFO, LDC, M, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION AP( * ), C( LDC, * ), TAU( * ), WORK( * ) 00015 * .. 00016 * 00017 * Purpose 00018 * ======= 00019 * 00020 * DOPMTR overwrites the general real M-by-N matrix C with 00021 * 00022 * SIDE = 'L' SIDE = 'R' 00023 * TRANS = 'N': Q * C C * Q 00024 * TRANS = 'T': Q**T * C C * Q**T 00025 * 00026 * where Q is a real orthogonal 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 DSPTRD 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**T from the Left; 00040 * = 'R': apply Q or Q**T from the Right. 00041 * 00042 * UPLO (input) CHARACTER*1 00043 * = 'U': Upper triangular packed storage used in previous 00044 * call to DSPTRD; 00045 * = 'L': Lower triangular packed storage used in previous 00046 * call to DSPTRD. 00047 * 00048 * TRANS (input) CHARACTER*1 00049 * = 'N': No transpose, apply Q; 00050 * = 'T': Transpose, apply Q**T. 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) DOUBLE PRECISION 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 DSPTRD. AP is modified by the routine but 00063 * restored on exit. 00064 * 00065 * TAU (input) DOUBLE PRECISION 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 DSPTRD. 00069 * 00070 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 00071 * On entry, the M-by-N matrix C. 00072 * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. 00073 * 00074 * LDC (input) INTEGER 00075 * The leading dimension of the array C. LDC >= max(1,M). 00076 * 00077 * WORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ONE 00089 PARAMETER ( ONE = 1.0D+0 ) 00090 * .. 00091 * .. Local Scalars .. 00092 LOGICAL FORWRD, LEFT, NOTRAN, UPPER 00093 INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ 00094 DOUBLE PRECISION AII 00095 * .. 00096 * .. External Functions .. 00097 LOGICAL LSAME 00098 EXTERNAL LSAME 00099 * .. 00100 * .. External Subroutines .. 00101 EXTERNAL DLARF, XERBLA 00102 * .. 00103 * .. Intrinsic Functions .. 00104 INTRINSIC 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, 'T' ) ) 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( 'DOPMTR', -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 DSPTRD 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) is applied to C(1:i,1:n) 00174 * 00175 MI = I 00176 ELSE 00177 * 00178 * H(i) is applied to C(1:m,1:i) 00179 * 00180 NI = I 00181 END IF 00182 * 00183 * Apply H(i) 00184 * 00185 AII = AP( II ) 00186 AP( II ) = ONE 00187 CALL DLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, LDC, 00188 $ WORK ) 00189 AP( II ) = AII 00190 * 00191 IF( FORWRD ) THEN 00192 II = II + I + 2 00193 ELSE 00194 II = II - I - 1 00195 END IF 00196 10 CONTINUE 00197 ELSE 00198 * 00199 * Q was determined by a call to DSPTRD with UPLO = 'L'. 00200 * 00201 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR. 00202 $ ( .NOT.LEFT .AND. NOTRAN ) 00203 * 00204 IF( FORWRD ) THEN 00205 I1 = 1 00206 I2 = NQ - 1 00207 I3 = 1 00208 II = 2 00209 ELSE 00210 I1 = NQ - 1 00211 I2 = 1 00212 I3 = -1 00213 II = NQ*( NQ+1 ) / 2 - 1 00214 END IF 00215 * 00216 IF( LEFT ) THEN 00217 NI = N 00218 JC = 1 00219 ELSE 00220 MI = M 00221 IC = 1 00222 END IF 00223 * 00224 DO 20 I = I1, I2, I3 00225 AII = AP( II ) 00226 AP( II ) = ONE 00227 IF( LEFT ) THEN 00228 * 00229 * H(i) is applied to C(i+1:m,1:n) 00230 * 00231 MI = M - I 00232 IC = I + 1 00233 ELSE 00234 * 00235 * H(i) is applied to C(1:m,i+1:n) 00236 * 00237 NI = N - I 00238 JC = I + 1 00239 END IF 00240 * 00241 * Apply H(i) 00242 * 00243 CALL DLARF( SIDE, MI, NI, AP( II ), 1, TAU( I ), 00244 $ C( IC, JC ), LDC, WORK ) 00245 AP( II ) = AII 00246 * 00247 IF( FORWRD ) THEN 00248 II = II + NQ - I + 1 00249 ELSE 00250 II = II - NQ + I - 2 00251 END IF 00252 20 CONTINUE 00253 END IF 00254 RETURN 00255 * 00256 * End of DOPMTR 00257 * 00258 END
1.7.3 
