Annotation of rpl/lapack/lapack/zupmtr.f, revision 1.3
1.1 bertrand 1: SUBROUTINE ZUPMTR( 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: 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)' is applied to C(1:i,1:n)
174: *
175: MI = I
176: ELSE
177: *
178: * H(i) or H(i)' is applied to C(1:m,1:i)
179: *
180: NI = I
181: END IF
182: *
183: * Apply H(i) or H(i)'
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)' 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)' 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)'
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
CVSweb interface <joel.bertrand@systella.fr>