LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
cunm22.f
Go to the documentation of this file.
1 *> \brief \b CUNM22 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 CUNM22 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cunm22.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cunm22.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cunm22.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CUNM22( 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 Q( LDQ, * ), C( LDC, * ), WORK( * )
30 * ..
31 *
32 *> \par Purpose
33 * ============
34 *>
35 *> \verbatim
36 *>
37 *> CUNM22 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 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 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 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 cunm22( 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 Q( ldq, * ), C( ldc, * ), WORK( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  COMPLEX ONE
184  parameter ( one = ( 1.0e+0, 0.0e+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 cgemm, clacpy, ctrmm, xerbla
196 * ..
197 * .. Intrinsic Functions ..
198  INTRINSIC cmplx, 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 ) = cmplx( lwkopt )
243  END IF
244 *
245  IF( info.NE.0 ) THEN
246  CALL xerbla( 'CUNM22', -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 CTRMM.
260 *
261  IF( n1.EQ.0 ) THEN
262  CALL ctrmm( 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 ctrmm( 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 clacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
286  $ ldwork )
287  CALL ctrmm( '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 cgemm( '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 clacpy( 'All', n2, len, c( 1, i ), ldc,
300  $ work( n1+1 ), ldwork )
301  CALL ctrmm( '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 cgemm( '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 clacpy( '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 clacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
324  $ ldwork )
325  CALL ctrmm( '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 cgemm( '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 clacpy( 'All', n1, len, c( 1, i ), ldc,
338  $ work( n2+1 ), ldwork )
339  CALL ctrmm( '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 cgemm( '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 clacpy( '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 clacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
364  $ ldwork )
365  CALL ctrmm( '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 cgemm( '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 clacpy( 'All', len, n1, c( i, 1 ), ldc,
378  $ work( 1 + n2*ldwork ), ldwork )
379  CALL ctrmm( '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 cgemm( '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 clacpy( '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 clacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
402  $ ldwork )
403  CALL ctrmm( '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 cgemm( '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 clacpy( 'All', len, n2, c( i, 1 ), ldc,
416  $ work( 1 + n1*ldwork ), ldwork )
417  CALL ctrmm( '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 cgemm( '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 clacpy( 'All', len, n, work, ldwork, c( i, 1 ),
430  $ ldc )
431  END DO
432  END IF
433  END IF
434 *
435  work( 1 ) = cmplx( lwkopt )
436  RETURN
437 *
438 * End of CUNM22
439 *
440  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:105
subroutine ctrmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRMM
Definition: ctrmm.f:179
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:189
subroutine cunm22(SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC, WORK, LWORK, INFO)
CUNM22 multiplies a general matrix by a banded unitary matrix.
Definition: cunm22.f:164