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
659