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