LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
cuncsd.f
Go to the documentation of this file.
1 *> \brief \b CUNCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CUNCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE CUNCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
22 * SIGNS, M, P, Q, X11, LDX11, X12,
23 * LDX12, X21, LDX21, X22, LDX22, THETA,
24 * U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
25 * LDV2T, WORK, LWORK, RWORK, LRWORK,
26 * IWORK, INFO )
27 *
28 * .. Scalar Arguments ..
29 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
30 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
31 * $ LDX21, LDX22, LRWORK, LWORK, M, P, Q
32 * ..
33 * .. Array Arguments ..
34 * INTEGER IWORK( * )
35 * REAL THETA( * )
36 * REAL RWORK( * )
37 * COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
38 * $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
39 * $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
40 * $ * )
41 * ..
42 *
43 *
44 *> \par Purpose:
45 * =============
46 *>
47 *> \verbatim
48 *>
49 *> CUNCSD computes the CS decomposition of an M-by-M partitioned
50 *> unitary matrix X:
51 *>
52 *> [ I 0 0 | 0 0 0 ]
53 *> [ 0 C 0 | 0 -S 0 ]
54 *> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**H
55 *> X = [-----------] = [---------] [---------------------] [---------] .
56 *> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ]
57 *> [ 0 S 0 | 0 C 0 ]
58 *> [ 0 0 I | 0 0 0 ]
59 *>
60 *> X11 is P-by-Q. The unitary matrices U1, U2, V1, and V2 are P-by-P,
61 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
62 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
63 *> which R = MIN(P,M-P,Q,M-Q).
64 *> \endverbatim
65 *
66 * Arguments:
67 * ==========
68 *
69 *> \param[in] JOBU1
70 *> \verbatim
71 *> JOBU1 is CHARACTER
72 *> = 'Y': U1 is computed;
73 *> otherwise: U1 is not computed.
74 *> \endverbatim
75 *>
76 *> \param[in] JOBU2
77 *> \verbatim
78 *> JOBU2 is CHARACTER
79 *> = 'Y': U2 is computed;
80 *> otherwise: U2 is not computed.
81 *> \endverbatim
82 *>
83 *> \param[in] JOBV1T
84 *> \verbatim
85 *> JOBV1T is CHARACTER
86 *> = 'Y': V1T is computed;
87 *> otherwise: V1T is not computed.
88 *> \endverbatim
89 *>
90 *> \param[in] JOBV2T
91 *> \verbatim
92 *> JOBV2T is CHARACTER
93 *> = 'Y': V2T is computed;
94 *> otherwise: V2T is not computed.
95 *> \endverbatim
96 *>
97 *> \param[in] TRANS
98 *> \verbatim
99 *> TRANS is CHARACTER
100 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
101 *> order;
102 *> otherwise: X, U1, U2, V1T, and V2T are stored in column-
103 *> major order.
104 *> \endverbatim
105 *>
106 *> \param[in] SIGNS
107 *> \verbatim
108 *> SIGNS is CHARACTER
109 *> = 'O': The lower-left block is made nonpositive (the
110 *> "other" convention);
111 *> otherwise: The upper-right block is made nonpositive (the
112 *> "default" convention).
113 *> \endverbatim
114 *>
115 *> \param[in] M
116 *> \verbatim
117 *> M is INTEGER
118 *> The number of rows and columns in X.
119 *> \endverbatim
120 *>
121 *> \param[in] P
122 *> \verbatim
123 *> P is INTEGER
124 *> The number of rows in X11 and X12. 0 <= P <= M.
125 *> \endverbatim
126 *>
127 *> \param[in] Q
128 *> \verbatim
129 *> Q is INTEGER
130 *> The number of columns in X11 and X21. 0 <= Q <= M.
131 *> \endverbatim
132 *>
133 *> \param[in,out] X11
134 *> \verbatim
135 *> X11 is COMPLEX array, dimension (LDX11,Q)
136 *> On entry, part of the unitary matrix whose CSD is desired.
137 *> \endverbatim
138 *>
139 *> \param[in] LDX11
140 *> \verbatim
141 *> LDX11 is INTEGER
142 *> The leading dimension of X11. LDX11 >= MAX(1,P).
143 *> \endverbatim
144 *>
145 *> \param[in,out] X12
146 *> \verbatim
147 *> X12 is COMPLEX array, dimension (LDX12,M-Q)
148 *> On entry, part of the unitary matrix whose CSD is desired.
149 *> \endverbatim
150 *>
151 *> \param[in] LDX12
152 *> \verbatim
153 *> LDX12 is INTEGER
154 *> The leading dimension of X12. LDX12 >= MAX(1,P).
155 *> \endverbatim
156 *>
157 *> \param[in,out] X21
158 *> \verbatim
159 *> X21 is COMPLEX array, dimension (LDX21,Q)
160 *> On entry, part of the unitary matrix whose CSD is desired.
161 *> \endverbatim
162 *>
163 *> \param[in] LDX21
164 *> \verbatim
165 *> LDX21 is INTEGER
166 *> The leading dimension of X11. LDX21 >= MAX(1,M-P).
167 *> \endverbatim
168 *>
169 *> \param[in,out] X22
170 *> \verbatim
171 *> X22 is COMPLEX array, dimension (LDX22,M-Q)
172 *> On entry, part of the unitary matrix whose CSD is desired.
173 *> \endverbatim
174 *>
175 *> \param[in] LDX22
176 *> \verbatim
177 *> LDX22 is INTEGER
178 *> The leading dimension of X11. LDX22 >= MAX(1,M-P).
179 *> \endverbatim
180 *>
181 *> \param[out] THETA
182 *> \verbatim
183 *> THETA is REAL array, dimension (R), in which R =
184 *> MIN(P,M-P,Q,M-Q).
185 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
186 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
187 *> \endverbatim
188 *>
189 *> \param[out] U1
190 *> \verbatim
191 *> U1 is COMPLEX array, dimension (LDU1,P)
192 *> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
193 *> \endverbatim
194 *>
195 *> \param[in] LDU1
196 *> \verbatim
197 *> LDU1 is INTEGER
198 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
199 *> MAX(1,P).
200 *> \endverbatim
201 *>
202 *> \param[out] U2
203 *> \verbatim
204 *> U2 is COMPLEX array, dimension (LDU2,M-P)
205 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
206 *> matrix U2.
207 *> \endverbatim
208 *>
209 *> \param[in] LDU2
210 *> \verbatim
211 *> LDU2 is INTEGER
212 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
213 *> MAX(1,M-P).
214 *> \endverbatim
215 *>
216 *> \param[out] V1T
217 *> \verbatim
218 *> V1T is COMPLEX array, dimension (LDV1T,Q)
219 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
220 *> matrix V1**H.
221 *> \endverbatim
222 *>
223 *> \param[in] LDV1T
224 *> \verbatim
225 *> LDV1T is INTEGER
226 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
227 *> MAX(1,Q).
228 *> \endverbatim
229 *>
230 *> \param[out] V2T
231 *> \verbatim
232 *> V2T is COMPLEX array, dimension (LDV2T,M-Q)
233 *> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) unitary
234 *> matrix V2**H.
235 *> \endverbatim
236 *>
237 *> \param[in] LDV2T
238 *> \verbatim
239 *> LDV2T is INTEGER
240 *> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
241 *> MAX(1,M-Q).
242 *> \endverbatim
243 *>
244 *> \param[out] WORK
245 *> \verbatim
246 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
247 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
248 *> \endverbatim
249 *>
250 *> \param[in] LWORK
251 *> \verbatim
252 *> LWORK is INTEGER
253 *> The dimension of the array WORK.
254 *>
255 *> If LWORK = -1, then a workspace query is assumed; the routine
256 *> only calculates the optimal size of the WORK array, returns
257 *> this value as the first entry of the work array, and no error
258 *> message related to LWORK is issued by XERBLA.
259 *> \endverbatim
260 *>
261 *> \param[out] RWORK
262 *> \verbatim
263 *> RWORK is REAL array, dimension MAX(1,LRWORK)
264 *> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
265 *> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
266 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
267 *> define the matrix in intermediate bidiagonal-block form
268 *> remaining after nonconvergence. INFO specifies the number
269 *> of nonzero PHI's.
270 *> \endverbatim
271 *>
272 *> \param[in] LRWORK
273 *> \verbatim
274 *> LRWORK is INTEGER
275 *> The dimension of the array RWORK.
276 *>
277 *> If LRWORK = -1, then a workspace query is assumed; the routine
278 *> only calculates the optimal size of the RWORK array, returns
279 *> this value as the first entry of the work array, and no error
280 *> message related to LRWORK is issued by XERBLA.
281 *> \endverbatim
282 *>
283 *> \param[out] IWORK
284 *> \verbatim
285 *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
286 *> \endverbatim
287 *>
288 *> \param[out] INFO
289 *> \verbatim
290 *> INFO is INTEGER
291 *> = 0: successful exit.
292 *> < 0: if INFO = -i, the i-th argument had an illegal value.
293 *> > 0: CBBCSD did not converge. See the description of RWORK
294 *> above for details.
295 *> \endverbatim
296 *
297 *> \par References:
298 * ================
299 *>
300 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
301 *> Algorithms, 50(1):33-65, 2009.
302 *
303 * Authors:
304 * ========
305 *
306 *> \author Univ. of Tennessee
307 *> \author Univ. of California Berkeley
308 *> \author Univ. of Colorado Denver
309 *> \author NAG Ltd.
310 *
311 *> \date June 2016
312 *
313 *> \ingroup complexOTHERcomputational
314 *
315 * =====================================================================
316  RECURSIVE SUBROUTINE cuncsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
317  $ SIGNS, M, P, Q, X11, LDX11, X12,
318  $ LDX12, X21, LDX21, X22, LDX22, THETA,
319  $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
320  $ LDV2T, WORK, LWORK, RWORK, LRWORK,
321  $ IWORK, INFO )
322 *
323 * -- LAPACK computational routine (version 3.7.1) --
324 * -- LAPACK is a software package provided by Univ. of Tennessee, --
325 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
326 * June 2016
327 *
328 * .. Scalar Arguments ..
329  CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
330  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
331  $ ldx21, ldx22, lrwork, lwork, m, p, q
332 * ..
333 * .. Array Arguments ..
334  INTEGER iwork( * )
335  REAL theta( * )
336  REAL rwork( * )
337  COMPLEX u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
338  $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
339  $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
340  $ * )
341 * ..
342 *
343 * ===================================================================
344 *
345 * .. Parameters ..
346  COMPLEX one, zero
347  PARAMETER ( one = (1.0e0,0.0e0),
348  $ zero = (0.0e0,0.0e0) )
349 * ..
350 * .. Local Scalars ..
351  CHARACTER transt, signst
352  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
353  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
354  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
355  $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
356  $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
357  $ lorbdbworkopt, lorglqwork, lorglqworkmin,
358  $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
359  $ lorgqrworkopt, lworkmin, lworkopt, p1, q1
360  LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
361  $ wantv1t, wantv2t
362  INTEGER lrworkmin, lrworkopt
363  LOGICAL lrquery
364 * ..
365 * .. External Subroutines ..
366  EXTERNAL xerbla, cbbcsd, clacpy, clapmr, clapmt,
367  $ cunbdb, cunglq, cungqr
368 * ..
369 * .. External Functions ..
370  LOGICAL lsame
371  EXTERNAL lsame
372 * ..
373 * .. Intrinsic Functions
374  INTRINSIC int, max, min
375 * ..
376 * .. Executable Statements ..
377 *
378 * Test input arguments
379 *
380  info = 0
381  wantu1 = lsame( jobu1, 'Y' )
382  wantu2 = lsame( jobu2, 'Y' )
383  wantv1t = lsame( jobv1t, 'Y' )
384  wantv2t = lsame( jobv2t, 'Y' )
385  colmajor = .NOT. lsame( trans, 'T' )
386  defaultsigns = .NOT. lsame( signs, 'O' )
387  lquery = lwork .EQ. -1
388  lrquery = lrwork .EQ. -1
389  IF( m .LT. 0 ) THEN
390  info = -7
391  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
392  info = -8
393  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
394  info = -9
395  ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
396  info = -11
397  ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
398  info = -11
399  ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
400  info = -13
401  ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
402  info = -13
403  ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
404  info = -15
405  ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
406  info = -15
407  ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
408  info = -17
409  ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
410  info = -17
411  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
412  info = -20
413  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
414  info = -22
415  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
416  info = -24
417  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
418  info = -26
419  END IF
420 *
421 * Work with transpose if convenient
422 *
423  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
424  IF( colmajor ) THEN
425  transt = 'T'
426  ELSE
427  transt = 'N'
428  END IF
429  IF( defaultsigns ) THEN
430  signst = 'O'
431  ELSE
432  signst = 'D'
433  END IF
434  CALL cuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
435  $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
436  $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
437  $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
438  $ info )
439  RETURN
440  END IF
441 *
442 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
443 * convenient
444 *
445  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
446  IF( defaultsigns ) THEN
447  signst = 'O'
448  ELSE
449  signst = 'D'
450  END IF
451  CALL cuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
452  $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
453  $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
454  $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
455  RETURN
456  END IF
457 *
458 * Compute workspace
459 *
460  IF( info .EQ. 0 ) THEN
461 *
462 * Real workspace
463 *
464  iphi = 2
465  ib11d = iphi + max( 1, q - 1 )
466  ib11e = ib11d + max( 1, q )
467  ib12d = ib11e + max( 1, q - 1 )
468  ib12e = ib12d + max( 1, q )
469  ib21d = ib12e + max( 1, q - 1 )
470  ib21e = ib21d + max( 1, q )
471  ib22d = ib21e + max( 1, q - 1 )
472  ib22e = ib22d + max( 1, q )
473  ibbcsd = ib22e + max( 1, q - 1 )
474  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
475  $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t,
476  $ v2t, ldv2t, theta, theta, theta, theta, theta,
477  $ theta, theta, theta, rwork, -1, childinfo )
478  lbbcsdworkopt = int( rwork(1) )
479  lbbcsdworkmin = lbbcsdworkopt
480  lrworkopt = ibbcsd + lbbcsdworkopt - 1
481  lrworkmin = ibbcsd + lbbcsdworkmin - 1
482  rwork(1) = lrworkopt
483 *
484 * Complex workspace
485 *
486  itaup1 = 2
487  itaup2 = itaup1 + max( 1, p )
488  itauq1 = itaup2 + max( 1, m - p )
489  itauq2 = itauq1 + max( 1, q )
490  iorgqr = itauq2 + max( 1, m - q )
491  CALL cungqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
492  $ childinfo )
493  lorgqrworkopt = int( work(1) )
494  lorgqrworkmin = max( 1, m - q )
495  iorglq = itauq2 + max( 1, m - q )
496  CALL cunglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
497  $ childinfo )
498  lorglqworkopt = int( work(1) )
499  lorglqworkmin = max( 1, m - q )
500  iorbdb = itauq2 + max( 1, m - q )
501  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
502  $ x21, ldx21, x22, ldx22, theta, theta, u1, u2,
503  $ v1t, v2t, work, -1, childinfo )
504  lorbdbworkopt = int( work(1) )
505  lorbdbworkmin = lorbdbworkopt
506  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
507  $ iorbdb + lorbdbworkopt ) - 1
508  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
509  $ iorbdb + lorbdbworkmin ) - 1
510  work(1) = max(lworkopt,lworkmin)
511 *
512  IF( lwork .LT. lworkmin
513  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
514  info = -22
515  ELSE IF( lrwork .LT. lrworkmin
516  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
517  info = -24
518  ELSE
519  lorgqrwork = lwork - iorgqr + 1
520  lorglqwork = lwork - iorglq + 1
521  lorbdbwork = lwork - iorbdb + 1
522  lbbcsdwork = lrwork - ibbcsd + 1
523  END IF
524  END IF
525 *
526 * Abort if any illegal arguments
527 *
528  IF( info .NE. 0 ) THEN
529  CALL xerbla( 'CUNCSD', -info )
530  RETURN
531  ELSE IF( lquery .OR. lrquery ) THEN
532  RETURN
533  END IF
534 *
535 * Transform to bidiagonal block form
536 *
537  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
538  $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
539  $ work(itaup2), work(itauq1), work(itauq2),
540  $ work(iorbdb), lorbdbwork, childinfo )
541 *
542 * Accumulate Householder reflectors
543 *
544  IF( colmajor ) THEN
545  IF( wantu1 .AND. p .GT. 0 ) THEN
546  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
547  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
548  $ lorgqrwork, info)
549  END IF
550  IF( wantu2 .AND. m-p .GT. 0 ) THEN
551  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
552  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
553  $ work(iorgqr), lorgqrwork, info )
554  END IF
555  IF( wantv1t .AND. q .GT. 0 ) THEN
556  CALL clacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
557  $ ldv1t )
558  v1t(1, 1) = one
559  DO j = 2, q
560  v1t(1,j) = zero
561  v1t(j,1) = zero
562  END DO
563  CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
564  $ work(iorglq), lorglqwork, info )
565  END IF
566  IF( wantv2t .AND. m-q .GT. 0 ) THEN
567  CALL clacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
568  IF( m-p .GT. q ) THEN
569  CALL clacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
570  $ v2t(p+1,p+1), ldv2t )
571  END IF
572  IF( m .GT. q ) THEN
573  CALL cunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
574  $ work(iorglq), lorglqwork, info )
575  END IF
576  END IF
577  ELSE
578  IF( wantu1 .AND. p .GT. 0 ) THEN
579  CALL clacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
580  CALL cunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
581  $ lorglqwork, info)
582  END IF
583  IF( wantu2 .AND. m-p .GT. 0 ) THEN
584  CALL clacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
585  CALL cunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
586  $ work(iorglq), lorglqwork, info )
587  END IF
588  IF( wantv1t .AND. q .GT. 0 ) THEN
589  CALL clacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
590  $ ldv1t )
591  v1t(1, 1) = one
592  DO j = 2, q
593  v1t(1,j) = zero
594  v1t(j,1) = zero
595  END DO
596  CALL cungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
597  $ work(iorgqr), lorgqrwork, info )
598  END IF
599  IF( wantv2t .AND. m-q .GT. 0 ) THEN
600  p1 = min( p+1, m )
601  q1 = min( q+1, m )
602  CALL clacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
603  IF ( m .GT. p+q ) THEN
604  CALL clacpy( 'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
605  $ v2t(p+1,p+1), ldv2t )
606  END IF
607  CALL cungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
608  $ work(iorgqr), lorgqrwork, info )
609  END IF
610  END IF
611 *
612 * Compute the CSD of the matrix in bidiagonal-block form
613 *
614  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
615  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
616  $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
617  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
618  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
619  $ lbbcsdwork, info )
620 *
621 * Permute rows and columns to place identity submatrices in top-
622 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
623 * block and/or bottom-right corner of (2,1)-block and/or top-left
624 * corner of (2,2)-block
625 *
626  IF( q .GT. 0 .AND. wantu2 ) THEN
627  DO i = 1, q
628  iwork(i) = m - p - q + i
629  END DO
630  DO i = q + 1, m - p
631  iwork(i) = i - q
632  END DO
633  IF( colmajor ) THEN
634  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
635  ELSE
636  CALL clapmr( .false., m-p, m-p, u2, ldu2, iwork )
637  END IF
638  END IF
639  IF( m .GT. 0 .AND. wantv2t ) THEN
640  DO i = 1, p
641  iwork(i) = m - p - q + i
642  END DO
643  DO i = p + 1, m - q
644  iwork(i) = i - p
645  END DO
646  IF( .NOT. colmajor ) THEN
647  CALL clapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
648  ELSE
649  CALL clapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
650  END IF
651  END IF
652 *
653  RETURN
654 *
655 * End CUNCSD
656 *
657  END
658 
cunglq
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
Definition: cunglq.f:129
clacpy
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
cuncsd
recursive subroutine cuncsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
CUNCSD
Definition: cuncsd.f:322
cungqr
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:130
clapmr
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: clapmr.f:106
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
clapmt
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: clapmt.f:106
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
cunbdb
subroutine cunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
CUNBDB
Definition: cunbdb.f:289
cbbcsd
subroutine cbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, RWORK, LRWORK, INFO)
CBBCSD
Definition: cbbcsd.f:334