SRC/cuncsd.f

   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 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.6.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, CLASCL,
 367      $                   CLASET, 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