1 SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER SIDE, TRANS, UPLO
11 INTEGER INFO, LDC, M, N
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 AP( * ), C( LDC, * ), TAU( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZUPMTR overwrites the general complex M-by-N matrix C with
21 *
22 * SIDE = 'L' SIDE = 'R'
23 * TRANS = 'N': Q * C C * Q
24 * TRANS = 'C': Q**H * C C * Q**H
25 *
26 * where Q is a complex unitary matrix of order nq, with nq = m if
27 * SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
28 * nq-1 elementary reflectors, as returned by ZHPTRD using packed
29 * storage:
30 *
31 * if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
32 *
33 * if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
34 *
35 * Arguments
36 * =========
37 *
38 * SIDE (input) CHARACTER*1
39 * = 'L': apply Q or Q**H from the Left;
40 * = 'R': apply Q or Q**H from the Right.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangular packed storage used in previous
44 * call to ZHPTRD;
45 * = 'L': Lower triangular packed storage used in previous
46 * call to ZHPTRD.
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': No transpose, apply Q;
50 * = 'C': Conjugate transpose, apply Q**H.
51 *
52 * M (input) INTEGER
53 * The number of rows of the matrix C. M >= 0.
54 *
55 * N (input) INTEGER
56 * The number of columns of the matrix C. N >= 0.
57 *
58 * AP (input) COMPLEX*16 array, dimension
59 * (M*(M+1)/2) if SIDE = 'L'
60 * (N*(N+1)/2) if SIDE = 'R'
61 * The vectors which define the elementary reflectors, as
62 * returned by ZHPTRD. AP is modified by the routine but
63 * restored on exit.
64 *
65 * TAU (input) COMPLEX*16 array, dimension (M-1) if SIDE = 'L'
66 * or (N-1) if SIDE = 'R'
67 * TAU(i) must contain the scalar factor of the elementary
68 * reflector H(i), as returned by ZHPTRD.
69 *
70 * C (input/output) COMPLEX*16 array, dimension (LDC,N)
71 * On entry, the M-by-N matrix C.
72 * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
73 *
74 * LDC (input) INTEGER
75 * The leading dimension of the array C. LDC >= max(1,M).
76 *
77 * WORK (workspace) COMPLEX*16 array, dimension
78 * (N) if SIDE = 'L'
79 * (M) if SIDE = 'R'
80 *
81 * INFO (output) INTEGER
82 * = 0: successful exit
83 * < 0: if INFO = -i, the i-th argument had an illegal value
84 *
85 * =====================================================================
86 *
87 * .. Parameters ..
88 COMPLEX*16 ONE
89 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL FORWRD, LEFT, NOTRAN, UPPER
93 INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ
94 COMPLEX*16 AII, TAUI
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL XERBLA, ZLARF
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC DCONJG, MAX
105 * ..
106 * .. Executable Statements ..
107 *
108 * Test the input arguments
109 *
110 INFO = 0
111 LEFT = LSAME( SIDE, 'L' )
112 NOTRAN = LSAME( TRANS, 'N' )
113 UPPER = LSAME( UPLO, 'U' )
114 *
115 * NQ is the order of Q
116 *
117 IF( LEFT ) THEN
118 NQ = M
119 ELSE
120 NQ = N
121 END IF
122 IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
123 INFO = -1
124 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
125 INFO = -2
126 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
127 INFO = -3
128 ELSE IF( M.LT.0 ) THEN
129 INFO = -4
130 ELSE IF( N.LT.0 ) THEN
131 INFO = -5
132 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
133 INFO = -9
134 END IF
135 IF( INFO.NE.0 ) THEN
136 CALL XERBLA( 'ZUPMTR', -INFO )
137 RETURN
138 END IF
139 *
140 * Quick return if possible
141 *
142 IF( M.EQ.0 .OR. N.EQ.0 )
143 $ RETURN
144 *
145 IF( UPPER ) THEN
146 *
147 * Q was determined by a call to ZHPTRD with UPLO = 'U'
148 *
149 FORWRD = ( LEFT .AND. NOTRAN ) .OR.
150 $ ( .NOT.LEFT .AND. .NOT.NOTRAN )
151 *
152 IF( FORWRD ) THEN
153 I1 = 1
154 I2 = NQ - 1
155 I3 = 1
156 II = 2
157 ELSE
158 I1 = NQ - 1
159 I2 = 1
160 I3 = -1
161 II = NQ*( NQ+1 ) / 2 - 1
162 END IF
163 *
164 IF( LEFT ) THEN
165 NI = N
166 ELSE
167 MI = M
168 END IF
169 *
170 DO 10 I = I1, I2, I3
171 IF( LEFT ) THEN
172 *
173 * H(i) or H(i)**H is applied to C(1:i,1:n)
174 *
175 MI = I
176 ELSE
177 *
178 * H(i) or H(i)**H is applied to C(1:m,1:i)
179 *
180 NI = I
181 END IF
182 *
183 * Apply H(i) or H(i)**H
184 *
185 IF( NOTRAN ) THEN
186 TAUI = TAU( I )
187 ELSE
188 TAUI = DCONJG( TAU( I ) )
189 END IF
190 AII = AP( II )
191 AP( II ) = ONE
192 CALL ZLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, LDC,
193 $ WORK )
194 AP( II ) = AII
195 *
196 IF( FORWRD ) THEN
197 II = II + I + 2
198 ELSE
199 II = II - I - 1
200 END IF
201 10 CONTINUE
202 ELSE
203 *
204 * Q was determined by a call to ZHPTRD with UPLO = 'L'.
205 *
206 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
207 $ ( .NOT.LEFT .AND. NOTRAN )
208 *
209 IF( FORWRD ) THEN
210 I1 = 1
211 I2 = NQ - 1
212 I3 = 1
213 II = 2
214 ELSE
215 I1 = NQ - 1
216 I2 = 1
217 I3 = -1
218 II = NQ*( NQ+1 ) / 2 - 1
219 END IF
220 *
221 IF( LEFT ) THEN
222 NI = N
223 JC = 1
224 ELSE
225 MI = M
226 IC = 1
227 END IF
228 *
229 DO 20 I = I1, I2, I3
230 AII = AP( II )
231 AP( II ) = ONE
232 IF( LEFT ) THEN
233 *
234 * H(i) or H(i)**H is applied to C(i+1:m,1:n)
235 *
236 MI = M - I
237 IC = I + 1
238 ELSE
239 *
240 * H(i) or H(i)**H is applied to C(1:m,i+1:n)
241 *
242 NI = N - I
243 JC = I + 1
244 END IF
245 *
246 * Apply H(i) or H(i)**H
247 *
248 IF( NOTRAN ) THEN
249 TAUI = TAU( I )
250 ELSE
251 TAUI = DCONJG( TAU( I ) )
252 END IF
253 CALL ZLARF( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, JC ),
254 $ LDC, WORK )
255 AP( II ) = AII
256 *
257 IF( FORWRD ) THEN
258 II = II + NQ - I + 1
259 ELSE
260 II = II - NQ + I - 2
261 END IF
262 20 CONTINUE
263 END IF
264 RETURN
265 *
266 * End of ZUPMTR
267 *
268 END
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER SIDE, TRANS, UPLO
11 INTEGER INFO, LDC, M, N
12 * ..
13 * .. Array Arguments ..
14 COMPLEX*16 AP( * ), C( LDC, * ), TAU( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * ZUPMTR overwrites the general complex M-by-N matrix C with
21 *
22 * SIDE = 'L' SIDE = 'R'
23 * TRANS = 'N': Q * C C * Q
24 * TRANS = 'C': Q**H * C C * Q**H
25 *
26 * where Q is a complex unitary matrix of order nq, with nq = m if
27 * SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
28 * nq-1 elementary reflectors, as returned by ZHPTRD using packed
29 * storage:
30 *
31 * if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
32 *
33 * if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
34 *
35 * Arguments
36 * =========
37 *
38 * SIDE (input) CHARACTER*1
39 * = 'L': apply Q or Q**H from the Left;
40 * = 'R': apply Q or Q**H from the Right.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangular packed storage used in previous
44 * call to ZHPTRD;
45 * = 'L': Lower triangular packed storage used in previous
46 * call to ZHPTRD.
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': No transpose, apply Q;
50 * = 'C': Conjugate transpose, apply Q**H.
51 *
52 * M (input) INTEGER
53 * The number of rows of the matrix C. M >= 0.
54 *
55 * N (input) INTEGER
56 * The number of columns of the matrix C. N >= 0.
57 *
58 * AP (input) COMPLEX*16 array, dimension
59 * (M*(M+1)/2) if SIDE = 'L'
60 * (N*(N+1)/2) if SIDE = 'R'
61 * The vectors which define the elementary reflectors, as
62 * returned by ZHPTRD. AP is modified by the routine but
63 * restored on exit.
64 *
65 * TAU (input) COMPLEX*16 array, dimension (M-1) if SIDE = 'L'
66 * or (N-1) if SIDE = 'R'
67 * TAU(i) must contain the scalar factor of the elementary
68 * reflector H(i), as returned by ZHPTRD.
69 *
70 * C (input/output) COMPLEX*16 array, dimension (LDC,N)
71 * On entry, the M-by-N matrix C.
72 * On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
73 *
74 * LDC (input) INTEGER
75 * The leading dimension of the array C. LDC >= max(1,M).
76 *
77 * WORK (workspace) COMPLEX*16 array, dimension
78 * (N) if SIDE = 'L'
79 * (M) if SIDE = 'R'
80 *
81 * INFO (output) INTEGER
82 * = 0: successful exit
83 * < 0: if INFO = -i, the i-th argument had an illegal value
84 *
85 * =====================================================================
86 *
87 * .. Parameters ..
88 COMPLEX*16 ONE
89 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
90 * ..
91 * .. Local Scalars ..
92 LOGICAL FORWRD, LEFT, NOTRAN, UPPER
93 INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ
94 COMPLEX*16 AII, TAUI
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL XERBLA, ZLARF
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC DCONJG, MAX
105 * ..
106 * .. Executable Statements ..
107 *
108 * Test the input arguments
109 *
110 INFO = 0
111 LEFT = LSAME( SIDE, 'L' )
112 NOTRAN = LSAME( TRANS, 'N' )
113 UPPER = LSAME( UPLO, 'U' )
114 *
115 * NQ is the order of Q
116 *
117 IF( LEFT ) THEN
118 NQ = M
119 ELSE
120 NQ = N
121 END IF
122 IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
123 INFO = -1
124 ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
125 INFO = -2
126 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
127 INFO = -3
128 ELSE IF( M.LT.0 ) THEN
129 INFO = -4
130 ELSE IF( N.LT.0 ) THEN
131 INFO = -5
132 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
133 INFO = -9
134 END IF
135 IF( INFO.NE.0 ) THEN
136 CALL XERBLA( 'ZUPMTR', -INFO )
137 RETURN
138 END IF
139 *
140 * Quick return if possible
141 *
142 IF( M.EQ.0 .OR. N.EQ.0 )
143 $ RETURN
144 *
145 IF( UPPER ) THEN
146 *
147 * Q was determined by a call to ZHPTRD with UPLO = 'U'
148 *
149 FORWRD = ( LEFT .AND. NOTRAN ) .OR.
150 $ ( .NOT.LEFT .AND. .NOT.NOTRAN )
151 *
152 IF( FORWRD ) THEN
153 I1 = 1
154 I2 = NQ - 1
155 I3 = 1
156 II = 2
157 ELSE
158 I1 = NQ - 1
159 I2 = 1
160 I3 = -1
161 II = NQ*( NQ+1 ) / 2 - 1
162 END IF
163 *
164 IF( LEFT ) THEN
165 NI = N
166 ELSE
167 MI = M
168 END IF
169 *
170 DO 10 I = I1, I2, I3
171 IF( LEFT ) THEN
172 *
173 * H(i) or H(i)**H is applied to C(1:i,1:n)
174 *
175 MI = I
176 ELSE
177 *
178 * H(i) or H(i)**H is applied to C(1:m,1:i)
179 *
180 NI = I
181 END IF
182 *
183 * Apply H(i) or H(i)**H
184 *
185 IF( NOTRAN ) THEN
186 TAUI = TAU( I )
187 ELSE
188 TAUI = DCONJG( TAU( I ) )
189 END IF
190 AII = AP( II )
191 AP( II ) = ONE
192 CALL ZLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, LDC,
193 $ WORK )
194 AP( II ) = AII
195 *
196 IF( FORWRD ) THEN
197 II = II + I + 2
198 ELSE
199 II = II - I - 1
200 END IF
201 10 CONTINUE
202 ELSE
203 *
204 * Q was determined by a call to ZHPTRD with UPLO = 'L'.
205 *
206 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
207 $ ( .NOT.LEFT .AND. NOTRAN )
208 *
209 IF( FORWRD ) THEN
210 I1 = 1
211 I2 = NQ - 1
212 I3 = 1
213 II = 2
214 ELSE
215 I1 = NQ - 1
216 I2 = 1
217 I3 = -1
218 II = NQ*( NQ+1 ) / 2 - 1
219 END IF
220 *
221 IF( LEFT ) THEN
222 NI = N
223 JC = 1
224 ELSE
225 MI = M
226 IC = 1
227 END IF
228 *
229 DO 20 I = I1, I2, I3
230 AII = AP( II )
231 AP( II ) = ONE
232 IF( LEFT ) THEN
233 *
234 * H(i) or H(i)**H is applied to C(i+1:m,1:n)
235 *
236 MI = M - I
237 IC = I + 1
238 ELSE
239 *
240 * H(i) or H(i)**H is applied to C(1:m,i+1:n)
241 *
242 NI = N - I
243 JC = I + 1
244 END IF
245 *
246 * Apply H(i) or H(i)**H
247 *
248 IF( NOTRAN ) THEN
249 TAUI = TAU( I )
250 ELSE
251 TAUI = DCONJG( TAU( I ) )
252 END IF
253 CALL ZLARF( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, JC ),
254 $ LDC, WORK )
255 AP( II ) = AII
256 *
257 IF( FORWRD ) THEN
258 II = II + NQ - I + 1
259 ELSE
260 II = II - NQ + I - 2
261 END IF
262 20 CONTINUE
263 END IF
264 RETURN
265 *
266 * End of ZUPMTR
267 *
268 END