LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zunm22.f
Go to the documentation of this file.
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
subroutine ztrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRMM
Definition: ztrmm.f:177
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
subroutine zunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
ZUNM22 multiplies a general matrix by a banded unitary matrix.
Definition: zunm22.f:162