1 SUBROUTINE DOPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK,
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER SIDE, TRANS, UPLO
11 INTEGER INFO, LDC, M, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AP( * ), C( LDC, * ), TAU( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DOPMTR overwrites the general real M-by-N matrix C with
21 *
22 * SIDE = 'L' SIDE = 'R'
23 * TRANS = 'N': Q * C C * Q
24 * TRANS = 'T': Q**T * C C * Q**T
25 *
26 * where Q is a real orthogonal 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 DSPTRD 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**T from the Left;
40 * = 'R': apply Q or Q**T from the Right.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangular packed storage used in previous
44 * call to DSPTRD;
45 * = 'L': Lower triangular packed storage used in previous
46 * call to DSPTRD.
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': No transpose, apply Q;
50 * = 'T': Transpose, apply Q**T.
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) DOUBLE PRECISION 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 DSPTRD. AP is modified by the routine but
63 * restored on exit.
64 *
65 * TAU (input) DOUBLE PRECISION 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 DSPTRD.
69 *
70 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
71 * On entry, the M-by-N matrix C.
72 * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
73 *
74 * LDC (input) INTEGER
75 * The leading dimension of the array C. LDC >= max(1,M).
76 *
77 * WORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ONE
89 PARAMETER ( ONE = 1.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 DOUBLE PRECISION AII
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL DLARF, XERBLA
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC 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, 'T' ) ) 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( 'DOPMTR', -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 DSPTRD 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) is applied to C(1:i,1:n)
174 *
175 MI = I
176 ELSE
177 *
178 * H(i) is applied to C(1:m,1:i)
179 *
180 NI = I
181 END IF
182 *
183 * Apply H(i)
184 *
185 AII = AP( II )
186 AP( II ) = ONE
187 CALL DLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, LDC,
188 $ WORK )
189 AP( II ) = AII
190 *
191 IF( FORWRD ) THEN
192 II = II + I + 2
193 ELSE
194 II = II - I - 1
195 END IF
196 10 CONTINUE
197 ELSE
198 *
199 * Q was determined by a call to DSPTRD with UPLO = 'L'.
200 *
201 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
202 $ ( .NOT.LEFT .AND. NOTRAN )
203 *
204 IF( FORWRD ) THEN
205 I1 = 1
206 I2 = NQ - 1
207 I3 = 1
208 II = 2
209 ELSE
210 I1 = NQ - 1
211 I2 = 1
212 I3 = -1
213 II = NQ*( NQ+1 ) / 2 - 1
214 END IF
215 *
216 IF( LEFT ) THEN
217 NI = N
218 JC = 1
219 ELSE
220 MI = M
221 IC = 1
222 END IF
223 *
224 DO 20 I = I1, I2, I3
225 AII = AP( II )
226 AP( II ) = ONE
227 IF( LEFT ) THEN
228 *
229 * H(i) is applied to C(i+1:m,1:n)
230 *
231 MI = M - I
232 IC = I + 1
233 ELSE
234 *
235 * H(i) is applied to C(1:m,i+1:n)
236 *
237 NI = N - I
238 JC = I + 1
239 END IF
240 *
241 * Apply H(i)
242 *
243 CALL DLARF( SIDE, MI, NI, AP( II ), 1, TAU( I ),
244 $ C( IC, JC ), LDC, WORK )
245 AP( II ) = AII
246 *
247 IF( FORWRD ) THEN
248 II = II + NQ - I + 1
249 ELSE
250 II = II - NQ + I - 2
251 END IF
252 20 CONTINUE
253 END IF
254 RETURN
255 *
256 * End of DOPMTR
257 *
258 END
2 $ INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * November 2006
8 *
9 * .. Scalar Arguments ..
10 CHARACTER SIDE, TRANS, UPLO
11 INTEGER INFO, LDC, M, N
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AP( * ), C( LDC, * ), TAU( * ), WORK( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DOPMTR overwrites the general real M-by-N matrix C with
21 *
22 * SIDE = 'L' SIDE = 'R'
23 * TRANS = 'N': Q * C C * Q
24 * TRANS = 'T': Q**T * C C * Q**T
25 *
26 * where Q is a real orthogonal 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 DSPTRD 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**T from the Left;
40 * = 'R': apply Q or Q**T from the Right.
41 *
42 * UPLO (input) CHARACTER*1
43 * = 'U': Upper triangular packed storage used in previous
44 * call to DSPTRD;
45 * = 'L': Lower triangular packed storage used in previous
46 * call to DSPTRD.
47 *
48 * TRANS (input) CHARACTER*1
49 * = 'N': No transpose, apply Q;
50 * = 'T': Transpose, apply Q**T.
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) DOUBLE PRECISION 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 DSPTRD. AP is modified by the routine but
63 * restored on exit.
64 *
65 * TAU (input) DOUBLE PRECISION 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 DSPTRD.
69 *
70 * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
71 * On entry, the M-by-N matrix C.
72 * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
73 *
74 * LDC (input) INTEGER
75 * The leading dimension of the array C. LDC >= max(1,M).
76 *
77 * WORK (workspace) DOUBLE PRECISION 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 DOUBLE PRECISION ONE
89 PARAMETER ( ONE = 1.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 DOUBLE PRECISION AII
95 * ..
96 * .. External Functions ..
97 LOGICAL LSAME
98 EXTERNAL LSAME
99 * ..
100 * .. External Subroutines ..
101 EXTERNAL DLARF, XERBLA
102 * ..
103 * .. Intrinsic Functions ..
104 INTRINSIC 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, 'T' ) ) 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( 'DOPMTR', -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 DSPTRD 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) is applied to C(1:i,1:n)
174 *
175 MI = I
176 ELSE
177 *
178 * H(i) is applied to C(1:m,1:i)
179 *
180 NI = I
181 END IF
182 *
183 * Apply H(i)
184 *
185 AII = AP( II )
186 AP( II ) = ONE
187 CALL DLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, LDC,
188 $ WORK )
189 AP( II ) = AII
190 *
191 IF( FORWRD ) THEN
192 II = II + I + 2
193 ELSE
194 II = II - I - 1
195 END IF
196 10 CONTINUE
197 ELSE
198 *
199 * Q was determined by a call to DSPTRD with UPLO = 'L'.
200 *
201 FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
202 $ ( .NOT.LEFT .AND. NOTRAN )
203 *
204 IF( FORWRD ) THEN
205 I1 = 1
206 I2 = NQ - 1
207 I3 = 1
208 II = 2
209 ELSE
210 I1 = NQ - 1
211 I2 = 1
212 I3 = -1
213 II = NQ*( NQ+1 ) / 2 - 1
214 END IF
215 *
216 IF( LEFT ) THEN
217 NI = N
218 JC = 1
219 ELSE
220 MI = M
221 IC = 1
222 END IF
223 *
224 DO 20 I = I1, I2, I3
225 AII = AP( II )
226 AP( II ) = ONE
227 IF( LEFT ) THEN
228 *
229 * H(i) is applied to C(i+1:m,1:n)
230 *
231 MI = M - I
232 IC = I + 1
233 ELSE
234 *
235 * H(i) is applied to C(1:m,i+1:n)
236 *
237 NI = N - I
238 JC = I + 1
239 END IF
240 *
241 * Apply H(i)
242 *
243 CALL DLARF( SIDE, MI, NI, AP( II ), 1, TAU( I ),
244 $ C( IC, JC ), LDC, WORK )
245 AP( II ) = AII
246 *
247 IF( FORWRD ) THEN
248 II = II + NQ - I + 1
249 ELSE
250 II = II - NQ + I - 2
251 END IF
252 20 CONTINUE
253 END IF
254 RETURN
255 *
256 * End of DOPMTR
257 *
258 END