1: *> \brief \b ZUNM22 multiplies a general matrix by a banded unitary matrix.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZUNM22 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zunm22.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zunm22.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zunm22.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZUNM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
22: * $ WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER SIDE, TRANS
26: * INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 Q( LDQ, * ), C( LDC, * ), WORK( * )
30: * ..
31: *
32: *> \par Purpose
33: * ============
34: *>
35: *> \verbatim
36: *>
37: *> ZUNM22 overwrites the general complex M-by-N matrix C with
38: *>
39: *> SIDE = 'L' SIDE = 'R'
40: *> TRANS = 'N': Q * C C * Q
41: *> TRANS = 'C': Q**H * C C * Q**H
42: *>
43: *> where Q is a complex unitary matrix of order NQ, with NQ = M if
44: *> SIDE = 'L' and NQ = N if SIDE = 'R'.
45: *> The unitary matrix Q processes a 2-by-2 block structure
46: *>
47: *> [ Q11 Q12 ]
48: *> Q = [ ]
49: *> [ Q21 Q22 ],
50: *>
51: *> where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
52: *> N2-by-N2 upper triangular matrix.
53: *> \endverbatim
54: *
55: * Arguments:
56: * ==========
57: *
58: *> \param[in] SIDE
59: *> \verbatim
60: *> SIDE is CHARACTER*1
61: *> = 'L': apply Q or Q**H from the Left;
62: *> = 'R': apply Q or Q**H from the Right.
63: *> \endverbatim
64: *>
65: *> \param[in] TRANS
66: *> \verbatim
67: *> TRANS is CHARACTER*1
68: *> = 'N': apply Q (No transpose);
69: *> = 'C': apply Q**H (Conjugate transpose).
70: *> \endverbatim
71: *>
72: *> \param[in] M
73: *> \verbatim
74: *> M is INTEGER
75: *> The number of rows of the matrix C. M >= 0.
76: *> \endverbatim
77: *>
78: *> \param[in] N
79: *> \verbatim
80: *> N is INTEGER
81: *> The number of columns of the matrix C. N >= 0.
82: *> \endverbatim
83: *>
84: *> \param[in] N1
85: *> \param[in] N2
86: *> \verbatim
87: *> N1 is INTEGER
88: *> N2 is INTEGER
89: *> The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
90: *> The following requirement must be satisfied:
91: *> N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
92: *> \endverbatim
93: *>
94: *> \param[in] Q
95: *> \verbatim
96: *> Q is COMPLEX*16 array, dimension
97: *> (LDQ,M) if SIDE = 'L'
98: *> (LDQ,N) if SIDE = 'R'
99: *> \endverbatim
100: *>
101: *> \param[in] LDQ
102: *> \verbatim
103: *> LDQ is INTEGER
104: *> The leading dimension of the array Q.
105: *> LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
106: *> \endverbatim
107: *>
108: *> \param[in,out] C
109: *> \verbatim
110: *> C is COMPLEX*16 array, dimension (LDC,N)
111: *> On entry, the M-by-N matrix C.
112: *> On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
113: *> \endverbatim
114: *>
115: *> \param[in] LDC
116: *> \verbatim
117: *> LDC is INTEGER
118: *> The leading dimension of the array C. LDC >= max(1,M).
119: *> \endverbatim
120: *>
121: *> \param[out] WORK
122: *> \verbatim
123: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
124: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
125: *> \endverbatim
126: *>
127: *> \param[in] LWORK
128: *> \verbatim
129: *> LWORK is INTEGER
130: *> The dimension of the array WORK.
131: *> If SIDE = 'L', LWORK >= max(1,N);
132: *> if SIDE = 'R', LWORK >= max(1,M).
133: *> For optimum performance LWORK >= M*N.
134: *>
135: *> If LWORK = -1, then a workspace query is assumed; the routine
136: *> only calculates the optimal size of the WORK array, returns
137: *> this value as the first entry of the WORK array, and no error
138: *> message related to LWORK is issued by XERBLA.
139: *> \endverbatim
140: *>
141: *> \param[out] INFO
142: *> \verbatim
143: *> INFO is INTEGER
144: *> = 0: successful exit
145: *> < 0: if INFO = -i, the i-th argument had an illegal value
146: *> \endverbatim
147: *
148: *
149: * Authors:
150: * ========
151: *
152: *> \author Univ. of Tennessee
153: *> \author Univ. of California Berkeley
154: *> \author Univ. of Colorado Denver
155: *> \author NAG Ltd.
156: *
157: *> \ingroup complexOTHERcomputational
158: *
159: * =====================================================================
160: SUBROUTINE ZUNM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
161: $ WORK, LWORK, INFO )
162: *
163: * -- LAPACK computational routine --
164: * -- LAPACK is a software package provided by Univ. of Tennessee, --
165: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166: *
167: IMPLICIT NONE
168: *
169: * .. Scalar Arguments ..
170: CHARACTER SIDE, TRANS
171: INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
172: * ..
173: * .. Array Arguments ..
174: COMPLEX*16 Q( LDQ, * ), C( LDC, * ), WORK( * )
175: * ..
176: *
177: * =====================================================================
178: *
179: * .. Parameters ..
180: COMPLEX*16 ONE
181: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
182: *
183: * .. Local Scalars ..
184: LOGICAL LEFT, LQUERY, NOTRAN
185: INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
186: * ..
187: * .. External Functions ..
188: LOGICAL LSAME
189: EXTERNAL LSAME
190: * ..
191: * .. External Subroutines ..
192: EXTERNAL ZGEMM, ZLACPY, ZTRMM, XERBLA
193: * ..
194: * .. Intrinsic Functions ..
195: INTRINSIC DCMPLX, MAX, MIN
196: * ..
197: * .. Executable Statements ..
198: *
199: * Test the input arguments
200: *
201: INFO = 0
202: LEFT = LSAME( SIDE, 'L' )
203: NOTRAN = LSAME( TRANS, 'N' )
204: LQUERY = ( LWORK.EQ.-1 )
205: *
206: * NQ is the order of Q;
207: * NW is the minimum dimension of WORK.
208: *
209: IF( LEFT ) THEN
210: NQ = M
211: ELSE
212: NQ = N
213: END IF
214: NW = NQ
215: IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
216: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
217: INFO = -1
218: ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
219: $ THEN
220: INFO = -2
221: ELSE IF( M.LT.0 ) THEN
222: INFO = -3
223: ELSE IF( N.LT.0 ) THEN
224: INFO = -4
225: ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
226: INFO = -5
227: ELSE IF( N2.LT.0 ) THEN
228: INFO = -6
229: ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
230: INFO = -8
231: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
232: INFO = -10
233: ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
234: INFO = -12
235: END IF
236: *
237: IF( INFO.EQ.0 ) THEN
238: LWKOPT = M*N
239: WORK( 1 ) = DCMPLX( LWKOPT )
240: END IF
241: *
242: IF( INFO.NE.0 ) THEN
243: CALL XERBLA( 'ZUNM22', -INFO )
244: RETURN
245: ELSE IF( LQUERY ) THEN
246: RETURN
247: END IF
248: *
249: * Quick return if possible
250: *
251: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
252: WORK( 1 ) = 1
253: RETURN
254: END IF
255: *
256: * Degenerate cases (N1 = 0 or N2 = 0) are handled using ZTRMM.
257: *
258: IF( N1.EQ.0 ) THEN
259: CALL ZTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
260: $ Q, LDQ, C, LDC )
261: WORK( 1 ) = ONE
262: RETURN
263: ELSE IF( N2.EQ.0 ) THEN
264: CALL ZTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
265: $ Q, LDQ, C, LDC )
266: WORK( 1 ) = ONE
267: RETURN
268: END IF
269: *
270: * Compute the largest chunk size available from the workspace.
271: *
272: NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
273: *
274: IF( LEFT ) THEN
275: IF( NOTRAN ) THEN
276: DO I = 1, N, NB
277: LEN = MIN( NB, N-I+1 )
278: LDWORK = M
279: *
280: * Multiply bottom part of C by Q12.
281: *
282: CALL ZLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
283: $ LDWORK )
284: CALL ZTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
285: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
286: $ LDWORK )
287: *
288: * Multiply top part of C by Q11.
289: *
290: CALL ZGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
291: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
292: $ LDWORK )
293: *
294: * Multiply top part of C by Q21.
295: *
296: CALL ZLACPY( 'All', N2, LEN, C( 1, I ), LDC,
297: $ WORK( N1+1 ), LDWORK )
298: CALL ZTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
299: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
300: $ WORK( N1+1 ), LDWORK )
301: *
302: * Multiply bottom part of C by Q22.
303: *
304: CALL ZGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
305: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
306: $ ONE, WORK( N1+1 ), LDWORK )
307: *
308: * Copy everything back.
309: *
310: CALL ZLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
311: $ LDC )
312: END DO
313: ELSE
314: DO I = 1, N, NB
315: LEN = MIN( NB, N-I+1 )
316: LDWORK = M
317: *
318: * Multiply bottom part of C by Q21**H.
319: *
320: CALL ZLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
321: $ LDWORK )
322: CALL ZTRMM( 'Left', 'Upper', 'Conjugate', 'Non-Unit',
323: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
324: $ LDWORK )
325: *
326: * Multiply top part of C by Q11**H.
327: *
328: CALL ZGEMM( 'Conjugate', 'No Transpose', N2, LEN, N1,
329: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
330: $ LDWORK )
331: *
332: * Multiply top part of C by Q12**H.
333: *
334: CALL ZLACPY( 'All', N1, LEN, C( 1, I ), LDC,
335: $ WORK( N2+1 ), LDWORK )
336: CALL ZTRMM( 'Left', 'Lower', 'Conjugate', 'Non-Unit',
337: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
338: $ WORK( N2+1 ), LDWORK )
339: *
340: * Multiply bottom part of C by Q22**H.
341: *
342: CALL ZGEMM( 'Conjugate', 'No Transpose', N1, LEN, N2,
343: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
344: $ ONE, WORK( N2+1 ), LDWORK )
345: *
346: * Copy everything back.
347: *
348: CALL ZLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
349: $ LDC )
350: END DO
351: END IF
352: ELSE
353: IF( NOTRAN ) THEN
354: DO I = 1, M, NB
355: LEN = MIN( NB, M-I+1 )
356: LDWORK = LEN
357: *
358: * Multiply right part of C by Q21.
359: *
360: CALL ZLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
361: $ LDWORK )
362: CALL ZTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
363: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
364: $ LDWORK )
365: *
366: * Multiply left part of C by Q11.
367: *
368: CALL ZGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
369: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
370: $ LDWORK )
371: *
372: * Multiply left part of C by Q12.
373: *
374: CALL ZLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
375: $ WORK( 1 + N2*LDWORK ), LDWORK )
376: CALL ZTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
377: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
378: $ WORK( 1 + N2*LDWORK ), LDWORK )
379: *
380: * Multiply right part of C by Q22.
381: *
382: CALL ZGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
383: $ ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
384: $ ONE, WORK( 1 + N2*LDWORK ), LDWORK )
385: *
386: * Copy everything back.
387: *
388: CALL ZLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
389: $ LDC )
390: END DO
391: ELSE
392: DO I = 1, M, NB
393: LEN = MIN( NB, M-I+1 )
394: LDWORK = LEN
395: *
396: * Multiply right part of C by Q12**H.
397: *
398: CALL ZLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
399: $ LDWORK )
400: CALL ZTRMM( 'Right', 'Lower', 'Conjugate', 'Non-Unit',
401: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
402: $ LDWORK )
403: *
404: * Multiply left part of C by Q11**H.
405: *
406: CALL ZGEMM( 'No Transpose', 'Conjugate', LEN, N1, N2,
407: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
408: $ LDWORK )
409: *
410: * Multiply left part of C by Q21**H.
411: *
412: CALL ZLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
413: $ WORK( 1 + N1*LDWORK ), LDWORK )
414: CALL ZTRMM( 'Right', 'Upper', 'Conjugate', 'Non-Unit',
415: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
416: $ WORK( 1 + N1*LDWORK ), LDWORK )
417: *
418: * Multiply right part of C by Q22**H.
419: *
420: CALL ZGEMM( 'No Transpose', 'Conjugate', LEN, N2, N1,
421: $ ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
422: $ ONE, WORK( 1 + N1*LDWORK ), LDWORK )
423: *
424: * Copy everything back.
425: *
426: CALL ZLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
427: $ LDC )
428: END DO
429: END IF
430: END IF
431: *
432: WORK( 1 ) = DCMPLX( LWKOPT )
433: RETURN
434: *
435: * End of ZUNM22
436: *
437: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zunm22.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack