LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 (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 (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 (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 (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 November 2013
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.5.0) --
324 * -- LAPACK is a software package provided by Univ. of Tennessee, --
325 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
326 * November 2013
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, clascl,
368 * ..
369 * .. External Functions ..
370  LOGICAL lsame
371  EXTERNAL lsame
372 * ..
373 * .. Intrinsic Functions
374  INTRINSIC cos, int, max, min, sin
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, 0, 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, 0, 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