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: *> \date January 2015
158: *
159: *> \ingroup complexOTHERcomputational
160: *
161: * =====================================================================
162: SUBROUTINE ZUNM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
163: $ WORK, LWORK, INFO )
164: *
165: * -- LAPACK computational routine (version 3.6.0) --
166: * -- LAPACK is a software package provided by Univ. of Tennessee, --
167: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168: * January 2015
169: *
170: IMPLICIT NONE
171: *
172: * .. Scalar Arguments ..
173: CHARACTER SIDE, TRANS
174: INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
175: * ..
176: * .. Array Arguments ..
177: COMPLEX*16 Q( LDQ, * ), C( LDC, * ), WORK( * )
178: * ..
179: *
180: * =====================================================================
181: *
182: * .. Parameters ..
183: COMPLEX*16 ONE
184: PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
185: *
186: * .. Local Scalars ..
187: LOGICAL LEFT, LQUERY, NOTRAN
188: INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
189: * ..
190: * .. External Functions ..
191: LOGICAL LSAME
192: EXTERNAL LSAME
193: * ..
194: * .. External Subroutines ..
195: EXTERNAL ZGEMM, ZLACPY, ZTRMM, XERBLA
196: * ..
197: * .. Intrinsic Functions ..
198: INTRINSIC DCMPLX, MAX, MIN
199: * ..
200: * .. Executable Statements ..
201: *
202: * Test the input arguments
203: *
204: INFO = 0
205: LEFT = LSAME( SIDE, 'L' )
206: NOTRAN = LSAME( TRANS, 'N' )
207: LQUERY = ( LWORK.EQ.-1 )
208: *
209: * NQ is the order of Q;
210: * NW is the minimum dimension of WORK.
211: *
212: IF( LEFT ) THEN
213: NQ = M
214: ELSE
215: NQ = N
216: END IF
217: NW = NQ
218: IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
219: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
220: INFO = -1
221: ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'C' ) )
222: $ THEN
223: INFO = -2
224: ELSE IF( M.LT.0 ) THEN
225: INFO = -3
226: ELSE IF( N.LT.0 ) THEN
227: INFO = -4
228: ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
229: INFO = -5
230: ELSE IF( N2.LT.0 ) THEN
231: INFO = -6
232: ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
233: INFO = -8
234: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
235: INFO = -10
236: ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
237: INFO = -12
238: END IF
239: *
240: IF( INFO.EQ.0 ) THEN
241: LWKOPT = M*N
242: WORK( 1 ) = DCMPLX( LWKOPT )
243: END IF
244: *
245: IF( INFO.NE.0 ) THEN
246: CALL XERBLA( 'ZUNM22', -INFO )
247: RETURN
248: ELSE IF( LQUERY ) THEN
249: RETURN
250: END IF
251: *
252: * Quick return if possible
253: *
254: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
255: WORK( 1 ) = 1
256: RETURN
257: END IF
258: *
259: * Degenerate cases (N1 = 0 or N2 = 0) are handled using ZTRMM.
260: *
261: IF( N1.EQ.0 ) THEN
262: CALL ZTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
263: $ Q, LDQ, C, LDC )
264: WORK( 1 ) = ONE
265: RETURN
266: ELSE IF( N2.EQ.0 ) THEN
267: CALL ZTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
268: $ Q, LDQ, C, LDC )
269: WORK( 1 ) = ONE
270: RETURN
271: END IF
272: *
273: * Compute the largest chunk size available from the workspace.
274: *
275: NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
276: *
277: IF( LEFT ) THEN
278: IF( NOTRAN ) THEN
279: DO I = 1, N, NB
280: LEN = MIN( NB, N-I+1 )
281: LDWORK = M
282: *
283: * Multiply bottom part of C by Q12.
284: *
285: CALL ZLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
286: $ LDWORK )
287: CALL ZTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
288: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
289: $ LDWORK )
290: *
291: * Multiply top part of C by Q11.
292: *
293: CALL ZGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
294: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
295: $ LDWORK )
296: *
297: * Multiply top part of C by Q21.
298: *
299: CALL ZLACPY( 'All', N2, LEN, C( 1, I ), LDC,
300: $ WORK( N1+1 ), LDWORK )
301: CALL ZTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
302: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
303: $ WORK( N1+1 ), LDWORK )
304: *
305: * Multiply bottom part of C by Q22.
306: *
307: CALL ZGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
308: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
309: $ ONE, WORK( N1+1 ), LDWORK )
310: *
311: * Copy everything back.
312: *
313: CALL ZLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
314: $ LDC )
315: END DO
316: ELSE
317: DO I = 1, N, NB
318: LEN = MIN( NB, N-I+1 )
319: LDWORK = M
320: *
321: * Multiply bottom part of C by Q21**H.
322: *
323: CALL ZLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
324: $ LDWORK )
325: CALL ZTRMM( 'Left', 'Upper', 'Conjugate', 'Non-Unit',
326: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
327: $ LDWORK )
328: *
329: * Multiply top part of C by Q11**H.
330: *
331: CALL ZGEMM( 'Conjugate', 'No Transpose', N2, LEN, N1,
332: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
333: $ LDWORK )
334: *
335: * Multiply top part of C by Q12**H.
336: *
337: CALL ZLACPY( 'All', N1, LEN, C( 1, I ), LDC,
338: $ WORK( N2+1 ), LDWORK )
339: CALL ZTRMM( 'Left', 'Lower', 'Conjugate', 'Non-Unit',
340: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
341: $ WORK( N2+1 ), LDWORK )
342: *
343: * Multiply bottom part of C by Q22**H.
344: *
345: CALL ZGEMM( 'Conjugate', 'No Transpose', N1, LEN, N2,
346: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
347: $ ONE, WORK( N2+1 ), LDWORK )
348: *
349: * Copy everything back.
350: *
351: CALL ZLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
352: $ LDC )
353: END DO
354: END IF
355: ELSE
356: IF( NOTRAN ) THEN
357: DO I = 1, M, NB
358: LEN = MIN( NB, M-I+1 )
359: LDWORK = LEN
360: *
361: * Multiply right part of C by Q21.
362: *
363: CALL ZLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
364: $ LDWORK )
365: CALL ZTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
366: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
367: $ LDWORK )
368: *
369: * Multiply left part of C by Q11.
370: *
371: CALL ZGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
372: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
373: $ LDWORK )
374: *
375: * Multiply left part of C by Q12.
376: *
377: CALL ZLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
378: $ WORK( 1 + N2*LDWORK ), LDWORK )
379: CALL ZTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
380: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
381: $ WORK( 1 + N2*LDWORK ), LDWORK )
382: *
383: * Multiply right part of C by Q22.
384: *
385: CALL ZGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
386: $ ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
387: $ ONE, WORK( 1 + N2*LDWORK ), LDWORK )
388: *
389: * Copy everything back.
390: *
391: CALL ZLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
392: $ LDC )
393: END DO
394: ELSE
395: DO I = 1, M, NB
396: LEN = MIN( NB, M-I+1 )
397: LDWORK = LEN
398: *
399: * Multiply right part of C by Q12**H.
400: *
401: CALL ZLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
402: $ LDWORK )
403: CALL ZTRMM( 'Right', 'Lower', 'Conjugate', 'Non-Unit',
404: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
405: $ LDWORK )
406: *
407: * Multiply left part of C by Q11**H.
408: *
409: CALL ZGEMM( 'No Transpose', 'Conjugate', LEN, N1, N2,
410: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
411: $ LDWORK )
412: *
413: * Multiply left part of C by Q21**H.
414: *
415: CALL ZLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
416: $ WORK( 1 + N1*LDWORK ), LDWORK )
417: CALL ZTRMM( 'Right', 'Upper', 'Conjugate', 'Non-Unit',
418: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
419: $ WORK( 1 + N1*LDWORK ), LDWORK )
420: *
421: * Multiply right part of C by Q22**H.
422: *
423: CALL ZGEMM( 'No Transpose', 'Conjugate', LEN, N2, N1,
424: $ ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
425: $ ONE, WORK( 1 + N1*LDWORK ), LDWORK )
426: *
427: * Copy everything back.
428: *
429: CALL ZLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
430: $ LDC )
431: END DO
432: END IF
433: END IF
434: *
435: WORK( 1 ) = DCMPLX( LWKOPT )
436: RETURN
437: *
438: * End of ZUNM22
439: *
440: 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