SRC/zuncsd2by1.f

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