LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 2011
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.4.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 2011
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
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) ) .OR.
396  $ ( .NOT.colmajor .AND. ldx11 .LT. max(1,q) ) ) THEN
397  info = -11
398  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
399  info = -20
400  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
401  info = -22
402  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
403  info = -24
404  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
405  info = -26
406  END IF
407 *
408 * Work with transpose if convenient
409 *
410  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
411  IF( colmajor ) THEN
412  transt = 'T'
413  ELSE
414  transt = 'N'
415  END IF
416  IF( defaultsigns ) THEN
417  signst = 'O'
418  ELSE
419  signst = 'D'
420  END IF
421  CALL cuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
422  $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
423  $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
424  $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
425  $ info )
426  return
427  END IF
428 *
429 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
430 * convenient
431 *
432  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
433  IF( defaultsigns ) THEN
434  signst = 'O'
435  ELSE
436  signst = 'D'
437  END IF
438  CALL cuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
439  $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
440  $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
441  $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
442  return
443  END IF
444 *
445 * Compute workspace
446 *
447  IF( info .EQ. 0 ) THEN
448 *
449 * Real workspace
450 *
451  iphi = 2
452  ib11d = iphi + max( 1, q - 1 )
453  ib11e = ib11d + max( 1, q )
454  ib12d = ib11e + max( 1, q - 1 )
455  ib12e = ib12d + max( 1, q )
456  ib21d = ib12e + max( 1, q - 1 )
457  ib21e = ib21d + max( 1, q )
458  ib22d = ib21e + max( 1, q - 1 )
459  ib22e = ib22d + max( 1, q )
460  ibbcsd = ib22e + max( 1, q - 1 )
461  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, 0,
462  $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, 0,
463  $ 0, 0, 0, 0, 0, 0, 0, rwork, -1, childinfo )
464  lbbcsdworkopt = int( rwork(1) )
465  lbbcsdworkmin = lbbcsdworkopt
466  lrworkopt = ibbcsd + lbbcsdworkopt - 1
467  lrworkmin = ibbcsd + lbbcsdworkmin - 1
468  rwork(1) = lrworkopt
469 *
470 * Complex workspace
471 *
472  itaup1 = 2
473  itaup2 = itaup1 + max( 1, p )
474  itauq1 = itaup2 + max( 1, m - p )
475  itauq2 = itauq1 + max( 1, q )
476  iorgqr = itauq2 + max( 1, m - q )
477  CALL cungqr( m-q, m-q, m-q, 0, max(1,m-q), 0, work, -1,
478  $ childinfo )
479  lorgqrworkopt = int( work(1) )
480  lorgqrworkmin = max( 1, m - q )
481  iorglq = itauq2 + max( 1, m - q )
482  CALL cunglq( m-q, m-q, m-q, 0, max(1,m-q), 0, work, -1,
483  $ childinfo )
484  lorglqworkopt = int( work(1) )
485  lorglqworkmin = max( 1, m - q )
486  iorbdb = itauq2 + max( 1, m - q )
487  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
488  $ x21, ldx21, x22, ldx22, 0, 0, 0, 0, 0, 0, work,
489  $ -1, childinfo )
490  lorbdbworkopt = int( work(1) )
491  lorbdbworkmin = lorbdbworkopt
492  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
493  $ iorbdb + lorbdbworkopt ) - 1
494  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
495  $ iorbdb + lorbdbworkmin ) - 1
496  work(1) = max(lworkopt,lworkmin)
497 *
498  IF( lwork .LT. lworkmin
499  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
500  info = -22
501  ELSE IF( lrwork .LT. lrworkmin
502  $ .AND. .NOT. ( lquery .OR. lrquery ) ) THEN
503  info = -24
504  ELSE
505  lorgqrwork = lwork - iorgqr + 1
506  lorglqwork = lwork - iorglq + 1
507  lorbdbwork = lwork - iorbdb + 1
508  lbbcsdwork = lrwork - ibbcsd + 1
509  END IF
510  END IF
511 *
512 * Abort if any illegal arguments
513 *
514  IF( info .NE. 0 ) THEN
515  CALL xerbla( 'CUNCSD', -info )
516  return
517  ELSE IF( lquery .OR. lrquery ) THEN
518  return
519  END IF
520 *
521 * Transform to bidiagonal block form
522 *
523  CALL cunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
524  $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
525  $ work(itaup2), work(itauq1), work(itauq2),
526  $ work(iorbdb), lorbdbwork, childinfo )
527 *
528 * Accumulate Householder reflectors
529 *
530  IF( colmajor ) THEN
531  IF( wantu1 .AND. p .GT. 0 ) THEN
532  CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
533  CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
534  $ lorgqrwork, info)
535  END IF
536  IF( wantu2 .AND. m-p .GT. 0 ) THEN
537  CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
538  CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
539  $ work(iorgqr), lorgqrwork, info )
540  END IF
541  IF( wantv1t .AND. q .GT. 0 ) THEN
542  CALL clacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
543  $ ldv1t )
544  v1t(1, 1) = one
545  DO j = 2, q
546  v1t(1,j) = zero
547  v1t(j,1) = zero
548  END DO
549  CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
550  $ work(iorglq), lorglqwork, info )
551  END IF
552  IF( wantv2t .AND. m-q .GT. 0 ) THEN
553  CALL clacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
554  CALL clacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
555  $ v2t(p+1,p+1), ldv2t )
556  CALL cunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
557  $ work(iorglq), lorglqwork, info )
558  END IF
559  ELSE
560  IF( wantu1 .AND. p .GT. 0 ) THEN
561  CALL clacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
562  CALL cunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
563  $ lorglqwork, info)
564  END IF
565  IF( wantu2 .AND. m-p .GT. 0 ) THEN
566  CALL clacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
567  CALL cunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
568  $ work(iorglq), lorglqwork, info )
569  END IF
570  IF( wantv1t .AND. q .GT. 0 ) THEN
571  CALL clacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
572  $ ldv1t )
573  v1t(1, 1) = one
574  DO j = 2, q
575  v1t(1,j) = zero
576  v1t(j,1) = zero
577  END DO
578  CALL cungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
579  $ work(iorgqr), lorgqrwork, info )
580  END IF
581  IF( wantv2t .AND. m-q .GT. 0 ) THEN
582  CALL clacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
583  CALL clacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
584  $ v2t(p+1,p+1), ldv2t )
585  CALL cungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
586  $ work(iorgqr), lorgqrwork, info )
587  END IF
588  END IF
589 *
590 * Compute the CSD of the matrix in bidiagonal-block form
591 *
592  CALL cbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
593  $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
594  $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
595  $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
596  $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
597  $ lbbcsdwork, info )
598 *
599 * Permute rows and columns to place identity submatrices in top-
600 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
601 * block and/or bottom-right corner of (2,1)-block and/or top-left
602 * corner of (2,2)-block
603 *
604  IF( q .GT. 0 .AND. wantu2 ) THEN
605  DO i = 1, q
606  iwork(i) = m - p - q + i
607  END DO
608  DO i = q + 1, m - p
609  iwork(i) = i - q
610  END DO
611  IF( colmajor ) THEN
612  CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
613  ELSE
614  CALL clapmr( .false., m-p, m-p, u2, ldu2, iwork )
615  END IF
616  END IF
617  IF( m .GT. 0 .AND. wantv2t ) THEN
618  DO i = 1, p
619  iwork(i) = m - p - q + i
620  END DO
621  DO i = p + 1, m - q
622  iwork(i) = i - p
623  END DO
624  IF( .NOT. colmajor ) THEN
625  CALL clapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
626  ELSE
627  CALL clapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
628  END IF
629  END IF
630 *
631  return
632 *
633 * End CUNCSD
634 *
635  END
636