SRC/cuncsd2by1.f

   1 *> \brief \b CUNCSD2BY1
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CUNCSD2BY1 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cuncsd2by1.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cuncsd2by1.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cuncsd2by1.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
  22 *                              X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
  23 *                              LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
  24 *                              INFO )
  25 *
  26 *       .. Scalar Arguments ..
  27 *       CHARACTER          JOBU1, JOBU2, JOBV1T
  28 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
  29 *      $                   M, P, Q
  30 *       INTEGER            LRWORK, LRWORKMIN, LRWORKOPT
  31 *       ..
  32 *       .. Array Arguments ..
  33 *       REAL               RWORK(*)
  34 *       REAL               THETA(*)
  35 *       COMPLEX            U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
  36 *      $                   X11(LDX11,*), X21(LDX21,*)
  37 *       INTEGER            IWORK(*)
  38 *       ..
  39 *
  40 *
  41 *> \par Purpose:
  42 *> =============
  43 *>
  44 *>\verbatim
  45 *>
  46 *> CUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
  47 *> orthonormal columns that has been partitioned into a 2-by-1 block
  48 *> structure:
  49 *>
  50 *>                                [  I1 0  0 ]
  51 *>                                [  0  C  0 ]
  52 *>          [ X11 ]   [ U1 |    ] [  0  0  0 ]
  53 *>      X = [-----] = [---------] [----------] V1**T .
  54 *>          [ X21 ]   [    | U2 ] [  0  0  0 ]
  55 *>                                [  0  S  0 ]
  56 *>                                [  0  0  I2]
  57 *>
  58 *> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P,
  59 *> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
  60 *> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
  61 *> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
  62 *> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
  63 *>
  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] M
  91 *> \verbatim
  92 *>          M is INTEGER
  93 *>          The number of rows in X.
  94 *> \endverbatim
  95 *>
  96 *> \param[in] P
  97 *> \verbatim
  98 *>          P is INTEGER
  99 *>          The number of rows in X11. 0 <= P <= M.
 100 *> \endverbatim
 101 *>
 102 *> \param[in] Q
 103 *> \verbatim
 104 *>          Q is INTEGER
 105 *>          The number of columns in X11 and X21. 0 <= Q <= M.
 106 *> \endverbatim
 107 *>
 108 *> \param[in,out] X11
 109 *> \verbatim
 110 *>          X11 is COMPLEX array, dimension (LDX11,Q)
 111 *>          On entry, part of the unitary matrix whose CSD is desired.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDX11
 115 *> \verbatim
 116 *>          LDX11 is INTEGER
 117 *>          The leading dimension of X11. LDX11 >= MAX(1,P).
 118 *> \endverbatim
 119 *>
 120 *> \param[in,out] X21
 121 *> \verbatim
 122 *>          X21 is COMPLEX array, dimension (LDX21,Q)
 123 *>          On entry, part of the unitary matrix whose CSD is desired.
 124 *> \endverbatim
 125 *>
 126 *> \param[in] LDX21
 127 *> \verbatim
 128 *>          LDX21 is INTEGER
 129 *>          The leading dimension of X21. LDX21 >= MAX(1,M-P).
 130 *> \endverbatim
 131 *>
 132 *> \param[out] THETA
 133 *> \verbatim
 134 *>          THETA is REAL array, dimension (R), in which R =
 135 *>          MIN(P,M-P,Q,M-Q).
 136 *>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
 137 *>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
 138 *> \endverbatim
 139 *>
 140 *> \param[out] U1
 141 *> \verbatim
 142 *>          U1 is COMPLEX array, dimension (P)
 143 *>          If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
 144 *> \endverbatim
 145 *>
 146 *> \param[in] LDU1
 147 *> \verbatim
 148 *>          LDU1 is INTEGER
 149 *>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
 150 *>          MAX(1,P).
 151 *> \endverbatim
 152 *>
 153 *> \param[out] U2
 154 *> \verbatim
 155 *>          U2 is COMPLEX array, dimension (M-P)
 156 *>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
 157 *>          matrix U2.
 158 *> \endverbatim
 159 *>
 160 *> \param[in] LDU2
 161 *> \verbatim
 162 *>          LDU2 is INTEGER
 163 *>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
 164 *>          MAX(1,M-P).
 165 *> \endverbatim
 166 *>
 167 *> \param[out] V1T
 168 *> \verbatim
 169 *>          V1T is COMPLEX array, dimension (Q)
 170 *>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
 171 *>          matrix V1**T.
 172 *> \endverbatim
 173 *>
 174 *> \param[in] LDV1T
 175 *> \verbatim
 176 *>          LDV1T is INTEGER
 177 *>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
 178 *>          MAX(1,Q).
 179 *> \endverbatim
 180 *>
 181 *> \param[out] WORK
 182 *> \verbatim
 183 *>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
 184 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 185 *> \endverbatim
 186 *>
 187 *> \param[in] LWORK
 188 *> \verbatim
 189 *>          LWORK is INTEGER
 190 *>          The dimension of the array WORK.
 191 *>
 192 *>          If LWORK = -1, then a workspace query is assumed; the routine
 193 *>          only calculates the optimal size of the WORK array, returns
 194 *>          this value as the first entry of the work array, and no error
 195 *>          message related to LWORK is issued by XERBLA.
 196 *> \endverbatim
 197 *>
 198 *> \param[out] RWORK
 199 *> \verbatim
 200 *>          RWORK is REAL array, dimension (MAX(1,LRWORK))
 201 *>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
 202 *>          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
 203 *>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
 204 *>          define the matrix in intermediate bidiagonal-block form
 205 *>          remaining after nonconvergence. INFO specifies the number
 206 *>          of nonzero PHI's.
 207 *> \endverbatim
 208 *>
 209 *> \param[in] LRWORK
 210 *> \verbatim
 211 *>          LRWORK is INTEGER
 212 *>          The dimension of the array RWORK.
 213 *>
 214 *>          If LRWORK = -1, then a workspace query is assumed; the routine
 215 *>          only calculates the optimal size of the RWORK array, returns
 216 *>          this value as the first entry of the work array, and no error
 217 *>          message related to LRWORK is issued by XERBLA.
 218 *> \endverbatim
 219 *
 220 *> \param[out] IWORK
 221 *> \verbatim
 222 *>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
 223 *> \endverbatim
 224 *>
 225 *> \param[out] INFO
 226 *> \verbatim
 227 *>          INFO is INTEGER
 228 *>          = 0:  successful exit.
 229 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 230 *>          > 0:  CBBCSD did not converge. See the description of WORK
 231 *>                above for details.
 232 *> \endverbatim
 233 *
 234 *> \par References:
 235 *  ================
 236 *>
 237 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
 238 *>      Algorithms, 50(1):33-65, 2009.
 239 *
 240 *  Authors:
 241 *  ========
 242 *
 243 *> \author Univ. of Tennessee
 244 *> \author Univ. of California Berkeley
 245 *> \author Univ. of Colorado Denver
 246 *> \author NAG Ltd.
 247 *
 248 *> \date June 2016
 249 *
 250 *> \ingroup complexOTHERcomputational
 251 *
 252 *  =====================================================================
 253       SUBROUTINE CUNCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
 254      $                       X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
 255      $                       LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
 256      $                       INFO )
 257 *
 258 *  -- LAPACK computational routine (version 3.6.1) --
 259 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 260 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 261 *     June 2016
 262 *
 263 *     .. Scalar Arguments ..
 264       CHARACTER          JOBU1, JOBU2, JOBV1T
 265       INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
 266      $                   M, P, Q
 267       INTEGER            LRWORK, LRWORKMIN, LRWORKOPT
 268 *     ..
 269 *     .. Array Arguments ..
 270       REAL               RWORK(*)
 271       REAL               THETA(*)
 272       COMPLEX            U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
 273      $                   X11(LDX11,*), X21(LDX21,*)
 274       INTEGER            IWORK(*)
 275 *     ..
 276 *
 277 *  =====================================================================
 278 *
 279 *     .. Parameters ..
 280       COMPLEX            ONE, ZERO
 281       PARAMETER          ( ONE = (1.0E0,0.0E0), ZERO = (0.0E0,0.0E0) )
 282 *     ..
 283 *     .. Local Scalars ..
 284       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
 285      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
 286      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
 287      $                   J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
 288      $                   LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
 289      $                   LWORKMIN, LWORKOPT, R
 290       LOGICAL            LQUERY, WANTU1, WANTU2, WANTV1T
 291 *     ..
 292 *     .. Local Arrays ..
 293       REAL               DUM( 1 )
 294       COMPLEX            CDUM( 1, 1 )
 295 *     ..
 296 *     .. External Subroutines ..
 297       EXTERNAL           CBBCSD, CCOPY, CLACPY, CLAPMR, CLAPMT, CUNBDB1,
 298      $                   CUNBDB2, CUNBDB3, CUNBDB4, CUNGLQ, CUNGQR,
 299      $                   XERBLA
 300 *     ..
 301 *     .. External Functions ..
 302       LOGICAL            LSAME
 303       EXTERNAL           LSAME
 304 *     ..
 305 *     .. Intrinsic Function ..
 306       INTRINSIC          INT, MAX, MIN
 307 *     ..
 308 *     .. Executable Statements ..
 309 *
 310 *     Test input arguments
 311 *
 312       INFO = 0
 313       WANTU1 = LSAME( JOBU1, 'Y' )
 314       WANTU2 = LSAME( JOBU2, 'Y' )
 315       WANTV1T = LSAME( JOBV1T, 'Y' )
 316       LQUERY = LWORK .EQ. -1
 317 *
 318       IF( M .LT. 0 ) THEN
 319          INFO = -4
 320       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
 321          INFO = -5
 322       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
 323          INFO = -6
 324       ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
 325          INFO = -8
 326       ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
 327          INFO = -10
 328       ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN
 329          INFO = -13
 330       ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN
 331          INFO = -15
 332       ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN
 333          INFO = -17
 334       END IF
 335 *
 336       R = MIN( P, M-P, Q, M-Q )
 337 *
 338 *     Compute workspace
 339 *
 340 *       WORK layout:
 341 *     |-----------------------------------------|
 342 *     | LWORKOPT (1)                            |
 343 *     |-----------------------------------------|
 344 *     | TAUP1 (MAX(1,P))                        |
 345 *     | TAUP2 (MAX(1,M-P))                      |
 346 *     | TAUQ1 (MAX(1,Q))                        |
 347 *     |-----------------------------------------|
 348 *     | CUNBDB WORK | CUNGQR WORK | CUNGLQ WORK |
 349 *     |             |             |             |
 350 *     |             |             |             |
 351 *     |             |             |             |
 352 *     |             |             |             |
 353 *     |-----------------------------------------|
 354 *       RWORK layout:
 355 *     |------------------|
 356 *     | LRWORKOPT (1)    |
 357 *     |------------------|
 358 *     | PHI (MAX(1,R-1)) |
 359 *     |------------------|
 360 *     | B11D (R)         |
 361 *     | B11E (R-1)       |
 362 *     | B12D (R)         |
 363 *     | B12E (R-1)       |
 364 *     | B21D (R)         |
 365 *     | B21E (R-1)       |
 366 *     | B22D (R)         |
 367 *     | B22E (R-1)       |
 368 *     | CBBCSD RWORK     |
 369 *     |------------------|
 370 *
 371       IF( INFO .EQ. 0 ) THEN
 372          IPHI = 2
 373          IB11D = IPHI + MAX( 1, R-1 )
 374          IB11E = IB11D + MAX( 1, R )
 375          IB12D = IB11E + MAX( 1, R - 1 )
 376          IB12E = IB12D + MAX( 1, R )
 377          IB21D = IB12E + MAX( 1, R - 1 )
 378          IB21E = IB21D + MAX( 1, R )
 379          IB22D = IB21E + MAX( 1, R - 1 )
 380          IB22E = IB22D + MAX( 1, R )
 381          IBBCSD = IB22E + MAX( 1, R - 1 )
 382          ITAUP1 = 2
 383          ITAUP2 = ITAUP1 + MAX( 1, P )
 384          ITAUQ1 = ITAUP2 + MAX( 1, M-P )
 385          IORBDB = ITAUQ1 + MAX( 1, Q )
 386          IORGQR = ITAUQ1 + MAX( 1, Q )
 387          IORGLQ = ITAUQ1 + MAX( 1, Q )
 388          LORGQRMIN = 1
 389          LORGQROPT = 1
 390          LORGLQMIN = 1
 391          LORGLQOPT = 1
 392          IF( R .EQ. Q ) THEN
 393             CALL CUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
 394      $                    DUM, CDUM, CDUM, CDUM, WORK, -1,
 395      $                    CHILDINFO )
 396             LORBDB = INT( WORK(1) )
 397             IF( WANTU1 .AND. P .GT. 0 ) THEN
 398                CALL CUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1,
 399      $                      CHILDINFO )
 400                LORGQRMIN = MAX( LORGQRMIN, P )
 401                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 402             ENDIF
 403             IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 404                CALL CUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1,
 405      $                      CHILDINFO )
 406                LORGQRMIN = MAX( LORGQRMIN, M-P )
 407                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 408             END IF
 409             IF( WANTV1T .AND. Q .GT. 0 ) THEN
 410                CALL CUNGLQ( Q-1, Q-1, Q-1, V1T, LDV1T,
 411      $                      CDUM, WORK(1), -1, CHILDINFO )
 412                LORGLQMIN = MAX( LORGLQMIN, Q-1 )
 413                LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
 414             END IF
 415             CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
 416      $                   DUM(1), U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM,
 417      $                   1, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
 418      $                   RWORK(1), -1, CHILDINFO )
 419             LBBCSD = INT( RWORK(1) )
 420          ELSE IF( R .EQ. P ) THEN
 421             CALL CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
 422      $                    CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO )
 423             LORBDB = INT( WORK(1) )
 424             IF( WANTU1 .AND. P .GT. 0 ) THEN
 425                CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, CDUM, WORK(1),
 426      $                      -1, CHILDINFO )
 427                LORGQRMIN = MAX( LORGQRMIN, P-1 )
 428                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 429             END IF
 430             IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 431                CALL CUNGQR( M-P, M-P, Q, U2, LDU2, CDUM, WORK(1), -1,
 432      $                      CHILDINFO )
 433                LORGQRMIN = MAX( LORGQRMIN, M-P )
 434                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 435             END IF
 436             IF( WANTV1T .AND. Q .GT. 0 ) THEN
 437                CALL CUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1,
 438      $                      CHILDINFO )
 439                LORGLQMIN = MAX( LORGLQMIN, Q )
 440                LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
 441             END IF
 442             CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
 443      $                   DUM, V1T, LDV1T, CDUM, 1, U1, LDU1, U2, LDU2,
 444      $                   DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
 445      $                   RWORK(1), -1, CHILDINFO )
 446             LBBCSD = INT( RWORK(1) )
 447          ELSE IF( R .EQ. M-P ) THEN
 448             CALL CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
 449      $                    CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO )
 450             LORBDB = INT( WORK(1) )
 451             IF( WANTU1 .AND. P .GT. 0 ) THEN
 452                CALL CUNGQR( P, P, Q, U1, LDU1, CDUM, WORK(1), -1,
 453      $                      CHILDINFO )
 454                LORGQRMIN = MAX( LORGQRMIN, P )
 455                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 456             END IF
 457             IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 458                CALL CUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2, CDUM,
 459      $                      WORK(1), -1, CHILDINFO )
 460                LORGQRMIN = MAX( LORGQRMIN, M-P-1 )
 461                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 462             END IF
 463             IF( WANTV1T .AND. Q .GT. 0 ) THEN
 464                CALL CUNGLQ( Q, Q, R, V1T, LDV1T, CDUM, WORK(1), -1,
 465      $                      CHILDINFO )
 466                LORGLQMIN = MAX( LORGLQMIN, Q )
 467                LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
 468             END IF
 469             CALL CBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
 470      $                   THETA, DUM, CDUM, 1, V1T, LDV1T, U2, LDU2, U1,
 471      $                   LDU1, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
 472      $                   RWORK(1), -1, CHILDINFO )
 473             LBBCSD = INT( RWORK(1) )
 474          ELSE
 475             CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DUM,
 476      $                    CDUM, CDUM, CDUM, CDUM, WORK(1), -1, CHILDINFO
 477      $                  )
 478             LORBDB = M + INT( WORK(1) )
 479             IF( WANTU1 .AND. P .GT. 0 ) THEN
 480                CALL CUNGQR( P, P, M-Q, U1, LDU1, CDUM, WORK(1), -1,
 481      $                      CHILDINFO )
 482                LORGQRMIN = MAX( LORGQRMIN, P )
 483                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 484             END IF
 485             IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 486                CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, CDUM, WORK(1), -1,
 487      $                      CHILDINFO )
 488                LORGQRMIN = MAX( LORGQRMIN, M-P )
 489                LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
 490             END IF
 491             IF( WANTV1T .AND. Q .GT. 0 ) THEN
 492                CALL CUNGLQ( Q, Q, Q, V1T, LDV1T, CDUM, WORK(1), -1,
 493      $                      CHILDINFO )
 494                LORGLQMIN = MAX( LORGLQMIN, Q )
 495                LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
 496             END IF
 497             CALL CBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
 498      $                   THETA, DUM, U2, LDU2, U1, LDU1, CDUM, 1, V1T,
 499      $                   LDV1T, DUM, DUM, DUM, DUM, DUM, DUM, DUM, DUM,
 500      $                   RWORK(1), -1, CHILDINFO )
 501             LBBCSD = INT( RWORK(1) )
 502          END IF
 503          LRWORKMIN = IBBCSD+LBBCSD-1
 504          LRWORKOPT = LRWORKMIN
 505          RWORK(1) = LRWORKOPT
 506          LWORKMIN = MAX( IORBDB+LORBDB-1,
 507      $                   IORGQR+LORGQRMIN-1,
 508      $                   IORGLQ+LORGLQMIN-1 )
 509          LWORKOPT = MAX( IORBDB+LORBDB-1,
 510      $                   IORGQR+LORGQROPT-1,
 511      $                   IORGLQ+LORGLQOPT-1 )
 512          WORK(1) = LWORKOPT
 513          IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
 514             INFO = -19
 515          END IF
 516       END IF
 517       IF( INFO .NE. 0 ) THEN
 518          CALL XERBLA( 'CUNCSD2BY1', -INFO )
 519          RETURN
 520       ELSE IF( LQUERY ) THEN
 521          RETURN
 522       END IF
 523       LORGQR = LWORK-IORGQR+1
 524       LORGLQ = LWORK-IORGLQ+1
 525 *
 526 *     Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
 527 *     in which R = MIN(P,M-P,Q,M-Q)
 528 *
 529       IF( R .EQ. Q ) THEN
 530 *
 531 *        Case 1: R = Q
 532 *
 533 *        Simultaneously bidiagonalize X11 and X21
 534 *
 535          CALL CUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
 536      $                 RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
 537      $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
 538 *
 539 *        Accumulate Householder reflectors
 540 *
 541          IF( WANTU1 .AND. P .GT. 0 ) THEN
 542             CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
 543             CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
 544      $                   LORGQR, CHILDINFO )
 545          END IF
 546          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 547             CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
 548             CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
 549      $                   WORK(IORGQR), LORGQR, CHILDINFO )
 550          END IF
 551          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 552             V1T(1,1) = ONE
 553             DO J = 2, Q
 554                V1T(1,J) = ZERO
 555                V1T(J,1) = ZERO
 556             END DO
 557             CALL CLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
 558      $                   LDV1T )
 559             CALL CUNGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
 560      $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
 561          END IF
 562 *
 563 *        Simultaneously diagonalize X11 and X21.
 564 *
 565          CALL CBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
 566      $                RWORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, CDUM,
 567      $                1, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D),
 568      $                RWORK(IB12E), RWORK(IB21D), RWORK(IB21E),
 569      $                RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD,
 570      $                CHILDINFO )
 571 *
 572 *        Permute rows and columns to place zero submatrices in
 573 *        preferred positions
 574 *
 575          IF( Q .GT. 0 .AND. WANTU2 ) THEN
 576             DO I = 1, Q
 577                IWORK(I) = M - P - Q + I
 578             END DO
 579             DO I = Q + 1, M - P
 580                IWORK(I) = I - Q
 581             END DO
 582             CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
 583          END IF
 584       ELSE IF( R .EQ. P ) THEN
 585 *
 586 *        Case 2: R = P
 587 *
 588 *        Simultaneously bidiagonalize X11 and X21
 589 *
 590          CALL CUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
 591      $                 RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
 592      $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
 593 *
 594 *        Accumulate Householder reflectors
 595 *
 596          IF( WANTU1 .AND. P .GT. 0 ) THEN
 597             U1(1,1) = ONE
 598             DO J = 2, P
 599                U1(1,J) = ZERO
 600                U1(J,1) = ZERO
 601             END DO
 602             CALL CLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
 603             CALL CUNGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
 604      $                   WORK(IORGQR), LORGQR, CHILDINFO )
 605          END IF
 606          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 607             CALL CLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
 608             CALL CUNGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
 609      $                   WORK(IORGQR), LORGQR, CHILDINFO )
 610          END IF
 611          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 612             CALL CLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T )
 613             CALL CUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
 614      $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
 615          END IF
 616 *
 617 *        Simultaneously diagonalize X11 and X21.
 618 *
 619          CALL CBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
 620      $                RWORK(IPHI), V1T, LDV1T, CDUM, 1, U1, LDU1, U2,
 621      $                LDU2, RWORK(IB11D), RWORK(IB11E), RWORK(IB12D),
 622      $                RWORK(IB12E), RWORK(IB21D), RWORK(IB21E),
 623      $                RWORK(IB22D), RWORK(IB22E), RWORK(IBBCSD), LBBCSD,
 624      $                CHILDINFO )
 625 *
 626 *        Permute rows and columns to place identity submatrices in
 627 *        preferred positions
 628 *
 629          IF( Q .GT. 0 .AND. WANTU2 ) THEN
 630             DO I = 1, Q
 631                IWORK(I) = M - P - Q + I
 632             END DO
 633             DO I = Q + 1, M - P
 634                IWORK(I) = I - Q
 635             END DO
 636             CALL CLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
 637          END IF
 638       ELSE IF( R .EQ. M-P ) THEN
 639 *
 640 *        Case 3: R = M-P
 641 *
 642 *        Simultaneously bidiagonalize X11 and X21
 643 *
 644          CALL CUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
 645      $                 RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
 646      $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
 647 *
 648 *        Accumulate Householder reflectors
 649 *
 650          IF( WANTU1 .AND. P .GT. 0 ) THEN
 651             CALL CLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
 652             CALL CUNGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
 653      $                   LORGQR, CHILDINFO )
 654          END IF
 655          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 656             U2(1,1) = ONE
 657             DO J = 2, M-P
 658                U2(1,J) = ZERO
 659                U2(J,1) = ZERO
 660             END DO
 661             CALL CLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
 662      $                   LDU2 )
 663             CALL CUNGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
 664      $                   WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
 665          END IF
 666          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 667             CALL CLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T )
 668             CALL CUNGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
 669      $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
 670          END IF
 671 *
 672 *        Simultaneously diagonalize X11 and X21.
 673 *
 674          CALL CBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
 675      $                THETA, RWORK(IPHI), CDUM, 1, V1T, LDV1T, U2, LDU2,
 676      $                U1, LDU1, RWORK(IB11D), RWORK(IB11E),
 677      $                RWORK(IB12D), RWORK(IB12E), RWORK(IB21D),
 678      $                RWORK(IB21E), RWORK(IB22D), RWORK(IB22E),
 679      $                RWORK(IBBCSD), LBBCSD, CHILDINFO )
 680 *
 681 *        Permute rows and columns to place identity submatrices in
 682 *        preferred positions
 683 *
 684          IF( Q .GT. R ) THEN
 685             DO I = 1, R
 686                IWORK(I) = Q - R + I
 687             END DO
 688             DO I = R + 1, Q
 689                IWORK(I) = I - R
 690             END DO
 691             IF( WANTU1 ) THEN
 692                CALL CLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
 693             END IF
 694             IF( WANTV1T ) THEN
 695                CALL CLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
 696             END IF
 697          END IF
 698       ELSE
 699 *
 700 *        Case 4: R = M-Q
 701 *
 702 *        Simultaneously bidiagonalize X11 and X21
 703 *
 704          CALL CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
 705      $                 RWORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
 706      $                 WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
 707      $                 LORBDB-M, CHILDINFO )
 708 *
 709 *        Accumulate Householder reflectors
 710 *
 711          IF( WANTU1 .AND. P .GT. 0 ) THEN
 712             CALL CCOPY( P, WORK(IORBDB), 1, U1, 1 )
 713             DO J = 2, P
 714                U1(1,J) = ZERO
 715             END DO
 716             CALL CLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
 717      $                   LDU1 )
 718             CALL CUNGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1),
 719      $                   WORK(IORGQR), LORGQR, CHILDINFO )
 720          END IF
 721          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 722             CALL CCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
 723             DO J = 2, M-P
 724                U2(1,J) = ZERO
 725             END DO
 726             CALL CLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
 727      $                   LDU2 )
 728             CALL CUNGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
 729      $                   WORK(IORGQR), LORGQR, CHILDINFO )
 730          END IF
 731          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 732             CALL CLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
 733             CALL CLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
 734      $                   V1T(M-Q+1,M-Q+1), LDV1T )
 735             CALL CLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
 736      $                   V1T(P+1,P+1), LDV1T )
 737             CALL CUNGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
 738      $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
 739          END IF
 740 *
 741 *        Simultaneously diagonalize X11 and X21.
 742 *
 743          CALL CBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
 744      $                THETA, RWORK(IPHI), U2, LDU2, U1, LDU1, CDUM, 1,
 745      $                V1T, LDV1T, RWORK(IB11D), RWORK(IB11E),
 746      $                RWORK(IB12D), RWORK(IB12E), RWORK(IB21D),
 747      $                RWORK(IB21E), RWORK(IB22D), RWORK(IB22E),
 748      $                RWORK(IBBCSD), LBBCSD, CHILDINFO )
 749 *
 750 *        Permute rows and columns to place identity submatrices in
 751 *        preferred positions
 752 *
 753          IF( P .GT. R ) THEN
 754             DO I = 1, R
 755                IWORK(I) = P - R + I
 756             END DO
 757             DO I = R + 1, P
 758                IWORK(I) = I - R
 759             END DO
 760             IF( WANTU1 ) THEN
 761                CALL CLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
 762             END IF
 763             IF( WANTV1T ) THEN
 764                CALL CLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
 765             END IF
 766          END IF
 767       END IF
 768 *
 769       RETURN
 770 *
 771 *     End of CUNCSD2BY1
 772 *
 773       END
 774