SRC/sorcsd.f

   1 *> \brief \b SORCSD
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SORCSD + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       RECURSIVE SUBROUTINE SORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
  22 *                                    SIGNS, M, P, Q, X11, LDX11, X12,
  23 *                                    LDX12, X21, LDX21, X22, LDX22, THETA,
  24 *                                    U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
  25 *                                    LDV2T, WORK, LWORK, IWORK, INFO )
  26 *
  27 *       .. Scalar Arguments ..
  28 *       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
  29 *       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
  30 *      $                   LDX21, LDX22, LWORK, M, P, Q
  31 *       ..
  32 *       .. Array Arguments ..
  33 *       INTEGER            IWORK( * )
  34 *       REAL               THETA( * )
  35 *       REAL               U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
  36 *      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
  37 *      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
  38 *      $                   * )
  39 *       ..
  40 *
  41 *
  42 *> \par Purpose:
  43 *  =============
  44 *>
  45 *> \verbatim
  46 *>
  47 *> SORCSD computes the CS decomposition of an M-by-M partitioned
  48 *> orthogonal matrix X:
  49 *>
  50 *>                                 [  I  0  0 |  0  0  0 ]
  51 *>                                 [  0  C  0 |  0 -S  0 ]
  52 *>     [ X11 | X12 ]   [ U1 |    ] [  0  0  0 |  0  0 -I ] [ V1 |    ]**T
  53 *> X = [-----------] = [---------] [---------------------] [---------]   .
  54 *>     [ X21 | X22 ]   [    | U2 ] [  0  0  0 |  I  0  0 ] [    | V2 ]
  55 *>                                 [  0  S  0 |  0  C  0 ]
  56 *>                                 [  0  0  I |  0  0  0 ]
  57 *>
  58 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
  59 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
  60 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
  61 *> which 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] JOBV2T
  89 *> \verbatim
  90 *>          JOBV2T is CHARACTER
  91 *>          = 'Y':      V2T is computed;
  92 *>          otherwise:  V2T is not computed.
  93 *> \endverbatim
  94 *>
  95 *> \param[in] TRANS
  96 *> \verbatim
  97 *>          TRANS is CHARACTER
  98 *>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
  99 *>                      order;
 100 *>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
 101 *>                      major order.
 102 *> \endverbatim
 103 *>
 104 *> \param[in] SIGNS
 105 *> \verbatim
 106 *>          SIGNS is CHARACTER
 107 *>          = 'O':      The lower-left block is made nonpositive (the
 108 *>                      "other" convention);
 109 *>          otherwise:  The upper-right block is made nonpositive (the
 110 *>                      "default" convention).
 111 *> \endverbatim
 112 *>
 113 *> \param[in] M
 114 *> \verbatim
 115 *>          M is INTEGER
 116 *>          The number of rows and columns in X.
 117 *> \endverbatim
 118 *>
 119 *> \param[in] P
 120 *> \verbatim
 121 *>          P is INTEGER
 122 *>          The number of rows in X11 and X12. 0 <= P <= M.
 123 *> \endverbatim
 124 *>
 125 *> \param[in] Q
 126 *> \verbatim
 127 *>          Q is INTEGER
 128 *>          The number of columns in X11 and X21. 0 <= Q <= M.
 129 *> \endverbatim
 130 *>
 131 *> \param[in,out] X11
 132 *> \verbatim
 133 *>          X11 is REAL array, dimension (LDX11,Q)
 134 *>          On entry, part of the orthogonal matrix whose CSD is desired.
 135 *> \endverbatim
 136 *>
 137 *> \param[in] LDX11
 138 *> \verbatim
 139 *>          LDX11 is INTEGER
 140 *>          The leading dimension of X11. LDX11 >= MAX(1,P).
 141 *> \endverbatim
 142 *>
 143 *> \param[in,out] X12
 144 *> \verbatim
 145 *>          X12 is REAL array, dimension (LDX12,M-Q)
 146 *>          On entry, part of the orthogonal matrix whose CSD is desired.
 147 *> \endverbatim
 148 *>
 149 *> \param[in] LDX12
 150 *> \verbatim
 151 *>          LDX12 is INTEGER
 152 *>          The leading dimension of X12. LDX12 >= MAX(1,P).
 153 *> \endverbatim
 154 *>
 155 *> \param[in,out] X21
 156 *> \verbatim
 157 *>          X21 is REAL array, dimension (LDX21,Q)
 158 *>          On entry, part of the orthogonal matrix whose CSD is desired.
 159 *> \endverbatim
 160 *>
 161 *> \param[in] LDX21
 162 *> \verbatim
 163 *>          LDX21 is INTEGER
 164 *>          The leading dimension of X11. LDX21 >= MAX(1,M-P).
 165 *> \endverbatim
 166 *>
 167 *> \param[in,out] X22
 168 *> \verbatim
 169 *>          X22 is REAL array, dimension (LDX22,M-Q)
 170 *>          On entry, part of the orthogonal matrix whose CSD is desired.
 171 *> \endverbatim
 172 *>
 173 *> \param[in] LDX22
 174 *> \verbatim
 175 *>          LDX22 is INTEGER
 176 *>          The leading dimension of X11. LDX22 >= MAX(1,M-P).
 177 *> \endverbatim
 178 *>
 179 *> \param[out] THETA
 180 *> \verbatim
 181 *>          THETA is REAL array, dimension (R), in which R =
 182 *>          MIN(P,M-P,Q,M-Q).
 183 *>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
 184 *>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
 185 *> \endverbatim
 186 *>
 187 *> \param[out] U1
 188 *> \verbatim
 189 *>          U1 is REAL array, dimension (P)
 190 *>          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
 191 *> \endverbatim
 192 *>
 193 *> \param[in] LDU1
 194 *> \verbatim
 195 *>          LDU1 is INTEGER
 196 *>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
 197 *>          MAX(1,P).
 198 *> \endverbatim
 199 *>
 200 *> \param[out] U2
 201 *> \verbatim
 202 *>          U2 is REAL array, dimension (M-P)
 203 *>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
 204 *>          matrix U2.
 205 *> \endverbatim
 206 *>
 207 *> \param[in] LDU2
 208 *> \verbatim
 209 *>          LDU2 is INTEGER
 210 *>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
 211 *>          MAX(1,M-P).
 212 *> \endverbatim
 213 *>
 214 *> \param[out] V1T
 215 *> \verbatim
 216 *>          V1T is REAL array, dimension (Q)
 217 *>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
 218 *>          matrix V1**T.
 219 *> \endverbatim
 220 *>
 221 *> \param[in] LDV1T
 222 *> \verbatim
 223 *>          LDV1T is INTEGER
 224 *>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
 225 *>          MAX(1,Q).
 226 *> \endverbatim
 227 *>
 228 *> \param[out] V2T
 229 *> \verbatim
 230 *>          V2T is REAL array, dimension (M-Q)
 231 *>          If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
 232 *>          matrix V2**T.
 233 *> \endverbatim
 234 *>
 235 *> \param[in] LDV2T
 236 *> \verbatim
 237 *>          LDV2T is INTEGER
 238 *>          The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
 239 *>          MAX(1,M-Q).
 240 *> \endverbatim
 241 *>
 242 *> \param[out] WORK
 243 *> \verbatim
 244 *>          WORK is REAL array, dimension (MAX(1,LWORK))
 245 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 246 *>          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
 247 *>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
 248 *>          define the matrix in intermediate bidiagonal-block form
 249 *>          remaining after nonconvergence. INFO specifies the number
 250 *>          of nonzero PHI's.
 251 *> \endverbatim
 252 *>
 253 *> \param[in] LWORK
 254 *> \verbatim
 255 *>          LWORK is INTEGER
 256 *>          The dimension of the array WORK.
 257 *>
 258 *>          If LWORK = -1, then a workspace query is assumed; the routine
 259 *>          only calculates the optimal size of the WORK array, returns
 260 *>          this value as the first entry of the work array, and no error
 261 *>          message related to LWORK is issued by XERBLA.
 262 *> \endverbatim
 263 *>
 264 *> \param[out] IWORK
 265 *> \verbatim
 266 *>          IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q))
 267 *> \endverbatim
 268 *>
 269 *> \param[out] INFO
 270 *> \verbatim
 271 *>          INFO is INTEGER
 272 *>          = 0:  successful exit.
 273 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 274 *>          > 0:  SBBCSD did not converge. See the description of WORK
 275 *>                above for details.
 276 *> \endverbatim
 277 *
 278 *> \par References:
 279 *  ================
 280 *>
 281 *>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
 282 *>      Algorithms, 50(1):33-65, 2009.
 283 *
 284 *  Authors:
 285 *  ========
 286 *
 287 *> \author Univ. of Tennessee
 288 *> \author Univ. of California Berkeley
 289 *> \author Univ. of Colorado Denver
 290 *> \author NAG Ltd.
 291 *
 292 *> \date November 2013
 293 *
 294 *> \ingroup realOTHERcomputational
 295 *
 296 *  =====================================================================
 297       RECURSIVE SUBROUTINE SORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
 298      $                             SIGNS, M, P, Q, X11, LDX11, X12,
 299      $                             LDX12, X21, LDX21, X22, LDX22, THETA,
 300      $                             U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
 301      $                             LDV2T, WORK, LWORK, IWORK, INFO )
 302 *
 303 *  -- LAPACK computational routine (version 3.5.0) --
 304 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 305 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 306 *     November 2013
 307 *
 308 *     .. Scalar Arguments ..
 309       CHARACTER          JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
 310       INTEGER            INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
 311      $                   LDX21, LDX22, LWORK, M, P, Q
 312 *     ..
 313 *     .. Array Arguments ..
 314       INTEGER            IWORK( * )
 315       REAL               THETA( * )
 316       REAL               U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
 317      $                   V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
 318      $                   X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
 319      $                   * )
 320 *     ..
 321 *
 322 *  ===================================================================
 323 *
 324 *     .. Parameters ..
 325       REAL               ONE, ZERO
 326       PARAMETER          ( ONE = 1.0E+0,
 327      $                     ZERO = 0.0E+0 )
 328 *     ..
 329 *     .. Local Arrays ..
 330       REAL               DUMMY(1)
 331 *     ..
 332 *     .. Local Scalars ..
 333       CHARACTER          TRANST, SIGNST
 334       INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
 335      $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
 336      $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
 337      $                   ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN,
 338      $                   LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN,
 339      $                   LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN,
 340      $                   LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN,
 341      $                   LORGQRWORKOPT, LWORKMIN, LWORKOPT
 342       LOGICAL            COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
 343      $                   WANTV1T, WANTV2T
 344 *     ..
 345 *     .. External Subroutines ..
 346       EXTERNAL           SBBCSD, SLACPY, SLAPMR, SLAPMT, SLASCL, SLASET,
 347      $                   SORBDB, SORGLQ, SORGQR, XERBLA
 348 *     ..
 349 *     .. External Functions ..
 350       LOGICAL            LSAME
 351       EXTERNAL           LSAME
 352 *     ..
 353 *     .. Intrinsic Functions
 354       INTRINSIC          INT, MAX, MIN
 355 *     ..
 356 *     .. Executable Statements ..
 357 *
 358 *     Test input arguments
 359 *
 360       INFO = 0
 361       WANTU1 = LSAME( JOBU1, 'Y' )
 362       WANTU2 = LSAME( JOBU2, 'Y' )
 363       WANTV1T = LSAME( JOBV1T, 'Y' )
 364       WANTV2T = LSAME( JOBV2T, 'Y' )
 365       COLMAJOR = .NOT. LSAME( TRANS, 'T' )
 366       DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' )
 367       LQUERY = LWORK .EQ. -1
 368       IF( M .LT. 0 ) THEN
 369          INFO = -7
 370       ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
 371          INFO = -8
 372       ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
 373          INFO = -9
 374       ELSE IF ( COLMAJOR .AND.  LDX11 .LT. MAX( 1, P ) ) THEN
 375         INFO = -11
 376       ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN
 377         INFO = -11
 378       ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN
 379         INFO = -13
 380       ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN
 381         INFO = -13
 382       ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN
 383         INFO = -15
 384       ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN
 385         INFO = -15
 386       ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN
 387         INFO = -17
 388       ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN
 389         INFO = -17
 390       ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN
 391          INFO = -20
 392       ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN
 393          INFO = -22
 394       ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN
 395          INFO = -24
 396       ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN
 397          INFO = -26
 398       END IF
 399 *
 400 *     Work with transpose if convenient
 401 *
 402       IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN
 403          IF( COLMAJOR ) THEN
 404             TRANST = 'T'
 405          ELSE
 406             TRANST = 'N'
 407          END IF
 408          IF( DEFAULTSIGNS ) THEN
 409             SIGNST = 'O'
 410          ELSE
 411             SIGNST = 'D'
 412          END IF
 413          CALL SORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M,
 414      $                Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22,
 415      $                LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1,
 416      $                U2, LDU2, WORK, LWORK, IWORK, INFO )
 417          RETURN
 418       END IF
 419 *
 420 *     Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
 421 *     convenient
 422 *
 423       IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN
 424          IF( DEFAULTSIGNS ) THEN
 425             SIGNST = 'O'
 426          ELSE
 427             SIGNST = 'D'
 428          END IF
 429          CALL SORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M,
 430      $                M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11,
 431      $                LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T,
 432      $                LDV1T, WORK, LWORK, IWORK, INFO )
 433          RETURN
 434       END IF
 435 *
 436 *     Compute workspace
 437 *
 438       IF( INFO .EQ. 0 ) THEN
 439 *
 440          IPHI = 2
 441          ITAUP1 = IPHI + MAX( 1, Q - 1 )
 442          ITAUP2 = ITAUP1 + MAX( 1, P )
 443          ITAUQ1 = ITAUP2 + MAX( 1, M - P )
 444          ITAUQ2 = ITAUQ1 + MAX( 1, Q )
 445          IORGQR = ITAUQ2 + MAX( 1, M - Q )
 446          CALL SORGQR( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1,
 447      $                CHILDINFO )
 448          LORGQRWORKOPT = INT( WORK(1) )
 449          LORGQRWORKMIN = MAX( 1, M - Q )
 450          IORGLQ = ITAUQ2 + MAX( 1, M - Q )
 451          CALL SORGLQ( M-Q, M-Q, M-Q, DUMMY, MAX(1,M-Q), DUMMY, WORK, -1,
 452      $                CHILDINFO )
 453          LORGLQWORKOPT = INT( WORK(1) )
 454          LORGLQWORKMIN = MAX( 1, M - Q )
 455          IORBDB = ITAUQ2 + MAX( 1, M - Q )
 456          CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12,
 457      $        X21, LDX21, X22, LDX22, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY,
 458      $        DUMMY,WORK,-1,CHILDINFO )
 459          LORBDBWORKOPT = INT( WORK(1) )
 460          LORBDBWORKMIN = LORBDBWORKOPT
 461          IB11D = ITAUQ2 + MAX( 1, M - Q )
 462          IB11E = IB11D + MAX( 1, Q )
 463          IB12D = IB11E + MAX( 1, Q - 1 )
 464          IB12E = IB12D + MAX( 1, Q )
 465          IB21D = IB12E + MAX( 1, Q - 1 )
 466          IB21E = IB21D + MAX( 1, Q )
 467          IB22D = IB21E + MAX( 1, Q - 1 )
 468          IB22E = IB22D + MAX( 1, Q )
 469          IBBCSD = IB22E + MAX( 1, Q - 1 )
 470          CALL SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
 471      $                DUMMY, DUMMY, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
 472      $                LDV2T, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY, DUMMY,
 473      $                DUMMY, DUMMY, WORK, -1, CHILDINFO )
 474          LBBCSDWORKOPT = INT( WORK(1) )
 475          LBBCSDWORKMIN = LBBCSDWORKOPT
 476          LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT,
 477      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1
 478          LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN,
 479      $              IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1
 480          WORK(1) = MAX(LWORKOPT,LWORKMIN)
 481 *
 482          IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN
 483             INFO = -22
 484          ELSE
 485             LORGQRWORK = LWORK - IORGQR + 1
 486             LORGLQWORK = LWORK - IORGLQ + 1
 487             LORBDBWORK = LWORK - IORBDB + 1
 488             LBBCSDWORK = LWORK - IBBCSD + 1
 489          END IF
 490       END IF
 491 *
 492 *     Abort if any illegal arguments
 493 *
 494       IF( INFO .NE. 0 ) THEN
 495          CALL XERBLA( 'SORCSD', -INFO )
 496          RETURN
 497       ELSE IF( LQUERY ) THEN
 498          RETURN
 499       END IF
 500 *
 501 *     Transform to bidiagonal block form
 502 *
 503       CALL SORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21,
 504      $             LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1),
 505      $             WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2),
 506      $             WORK(IORBDB), LORBDBWORK, CHILDINFO )
 507 *
 508 *     Accumulate Householder reflectors
 509 *
 510       IF( COLMAJOR ) THEN
 511          IF( WANTU1 .AND. P .GT. 0 ) THEN
 512             CALL SLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
 513             CALL SORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
 514      $                   LORGQRWORK, INFO)
 515          END IF
 516          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 517             CALL SLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
 518             CALL SORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
 519      $                   WORK(IORGQR), LORGQRWORK, INFO )
 520          END IF
 521          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 522             CALL SLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2),
 523      $                   LDV1T )
 524             V1T(1, 1) = ONE
 525             DO J = 2, Q
 526                V1T(1,J) = ZERO
 527                V1T(J,1) = ZERO
 528             END DO
 529             CALL SORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
 530      $                   WORK(IORGLQ), LORGLQWORK, INFO )
 531          END IF
 532          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
 533             CALL SLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T )
 534             CALL SLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22,
 535      $                   V2T(P+1,P+1), LDV2T )
 536             CALL SORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
 537      $                   WORK(IORGLQ), LORGLQWORK, INFO )
 538          END IF
 539       ELSE
 540          IF( WANTU1 .AND. P .GT. 0 ) THEN
 541             CALL SLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 )
 542             CALL SORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ),
 543      $                   LORGLQWORK, INFO)
 544          END IF
 545          IF( WANTU2 .AND. M-P .GT. 0 ) THEN
 546             CALL SLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 )
 547             CALL SORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
 548      $                   WORK(IORGLQ), LORGLQWORK, INFO )
 549          END IF
 550          IF( WANTV1T .AND. Q .GT. 0 ) THEN
 551             CALL SLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2),
 552      $                   LDV1T )
 553             V1T(1, 1) = ONE
 554             DO J = 2, Q
 555                V1T(1,J) = ZERO
 556                V1T(J,1) = ZERO
 557             END DO
 558             CALL SORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
 559      $                   WORK(IORGQR), LORGQRWORK, INFO )
 560          END IF
 561          IF( WANTV2T .AND. M-Q .GT. 0 ) THEN
 562             CALL SLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T )
 563             CALL SLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22,
 564      $                   V2T(P+1,P+1), LDV2T )
 565             CALL SORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2),
 566      $                   WORK(IORGQR), LORGQRWORK, INFO )
 567          END IF
 568       END IF
 569 *
 570 *     Compute the CSD of the matrix in bidiagonal-block form
 571 *
 572       CALL SBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA,
 573      $             WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
 574      $             LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D),
 575      $             WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D),
 576      $             WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO )
 577 *
 578 *     Permute rows and columns to place identity submatrices in top-
 579 *     left corner of (1,1)-block and/or bottom-right corner of (1,2)-
 580 *     block and/or bottom-right corner of (2,1)-block and/or top-left
 581 *     corner of (2,2)-block
 582 *
 583       IF( Q .GT. 0 .AND. WANTU2 ) THEN
 584          DO I = 1, Q
 585             IWORK(I) = M - P - Q + I
 586          END DO
 587          DO I = Q + 1, M - P
 588             IWORK(I) = I - Q
 589          END DO
 590          IF( COLMAJOR ) THEN
 591             CALL SLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
 592          ELSE
 593             CALL SLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK )
 594          END IF
 595       END IF
 596       IF( M .GT. 0 .AND. WANTV2T ) THEN
 597          DO I = 1, P
 598             IWORK(I) = M - P - Q + I
 599          END DO
 600          DO I = P + 1, M - Q
 601             IWORK(I) = I - P
 602          END DO
 603          IF( .NOT. COLMAJOR ) THEN
 604             CALL SLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
 605          ELSE
 606             CALL SLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK )
 607          END IF
 608       END IF
 609 *
 610       RETURN
 611 *
 612 *     End SORCSD
 613 *
 614       END
 615