1*> \brief \b DORCSD 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download DORCSD + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* RECURSIVE SUBROUTINE DORCSD( 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* DOUBLE PRECISION THETA( * ) 35* DOUBLE PRECISION 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*> DORCSD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDU1,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 DOUBLE PRECISION array, dimension (LDU2,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 DOUBLE PRECISION array, dimension (LDV1T,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 DOUBLE PRECISION array, dimension (LDV2T,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 DOUBLE PRECISION 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: DBBCSD 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*> \ingroup doubleOTHERcomputational 293* 294* ===================================================================== 295 RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, 296 $ SIGNS, M, P, Q, X11, LDX11, X12, 297 $ LDX12, X21, LDX21, X22, LDX22, THETA, 298 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 299 $ LDV2T, WORK, LWORK, IWORK, INFO ) 300* 301* -- LAPACK computational routine -- 302* -- LAPACK is a software package provided by Univ. of Tennessee, -- 303* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 304* 305* .. Scalar Arguments .. 306 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS 307 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12, 308 $ LDX21, LDX22, LWORK, M, P, Q 309* .. 310* .. Array Arguments .. 311 INTEGER IWORK( * ) 312 DOUBLE PRECISION THETA( * ) 313 DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ), 314 $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ), 315 $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22, 316 $ * ) 317* .. 318* 319* =================================================================== 320* 321* .. Parameters .. 322 DOUBLE PRECISION ONE, ZERO 323 PARAMETER ( ONE = 1.0D0, 324 $ ZERO = 0.0D0 ) 325* .. 326* .. Local Scalars .. 327 CHARACTER TRANST, SIGNST 328 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E, 329 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB, 330 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1, 331 $ ITAUQ2, J, LBBCSDWORK, LBBCSDWORKMIN, 332 $ LBBCSDWORKOPT, LORBDBWORK, LORBDBWORKMIN, 333 $ LORBDBWORKOPT, LORGLQWORK, LORGLQWORKMIN, 334 $ LORGLQWORKOPT, LORGQRWORK, LORGQRWORKMIN, 335 $ LORGQRWORKOPT, LWORKMIN, LWORKOPT 336 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2, 337 $ WANTV1T, WANTV2T 338* .. 339* .. External Subroutines .. 340 EXTERNAL DBBCSD, DLACPY, DLAPMR, DLAPMT, 341 $ DORBDB, DORGLQ, DORGQR, XERBLA 342* .. 343* .. External Functions .. 344 LOGICAL LSAME 345 EXTERNAL LSAME 346* .. 347* .. Intrinsic Functions 348 INTRINSIC INT, MAX, MIN 349* .. 350* .. Executable Statements .. 351* 352* Test input arguments 353* 354 INFO = 0 355 WANTU1 = LSAME( JOBU1, 'Y' ) 356 WANTU2 = LSAME( JOBU2, 'Y' ) 357 WANTV1T = LSAME( JOBV1T, 'Y' ) 358 WANTV2T = LSAME( JOBV2T, 'Y' ) 359 COLMAJOR = .NOT. LSAME( TRANS, 'T' ) 360 DEFAULTSIGNS = .NOT. LSAME( SIGNS, 'O' ) 361 LQUERY = LWORK .EQ. -1 362 IF( M .LT. 0 ) THEN 363 INFO = -7 364 ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN 365 INFO = -8 366 ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN 367 INFO = -9 368 ELSE IF ( COLMAJOR .AND. LDX11 .LT. MAX( 1, P ) ) THEN 369 INFO = -11 370 ELSE IF (.NOT. COLMAJOR .AND. LDX11 .LT. MAX( 1, Q ) ) THEN 371 INFO = -11 372 ELSE IF (COLMAJOR .AND. LDX12 .LT. MAX( 1, P ) ) THEN 373 INFO = -13 374 ELSE IF (.NOT. COLMAJOR .AND. LDX12 .LT. MAX( 1, M-Q ) ) THEN 375 INFO = -13 376 ELSE IF (COLMAJOR .AND. LDX21 .LT. MAX( 1, M-P ) ) THEN 377 INFO = -15 378 ELSE IF (.NOT. COLMAJOR .AND. LDX21 .LT. MAX( 1, Q ) ) THEN 379 INFO = -15 380 ELSE IF (COLMAJOR .AND. LDX22 .LT. MAX( 1, M-P ) ) THEN 381 INFO = -17 382 ELSE IF (.NOT. COLMAJOR .AND. LDX22 .LT. MAX( 1, M-Q ) ) THEN 383 INFO = -17 384 ELSE IF( WANTU1 .AND. LDU1 .LT. P ) THEN 385 INFO = -20 386 ELSE IF( WANTU2 .AND. LDU2 .LT. M-P ) THEN 387 INFO = -22 388 ELSE IF( WANTV1T .AND. LDV1T .LT. Q ) THEN 389 INFO = -24 390 ELSE IF( WANTV2T .AND. LDV2T .LT. M-Q ) THEN 391 INFO = -26 392 END IF 393* 394* Work with transpose if convenient 395* 396 IF( INFO .EQ. 0 .AND. MIN( P, M-P ) .LT. MIN( Q, M-Q ) ) THEN 397 IF( COLMAJOR ) THEN 398 TRANST = 'T' 399 ELSE 400 TRANST = 'N' 401 END IF 402 IF( DEFAULTSIGNS ) THEN 403 SIGNST = 'O' 404 ELSE 405 SIGNST = 'D' 406 END IF 407 CALL DORCSD( JOBV1T, JOBV2T, JOBU1, JOBU2, TRANST, SIGNST, M, 408 $ Q, P, X11, LDX11, X21, LDX21, X12, LDX12, X22, 409 $ LDX22, THETA, V1T, LDV1T, V2T, LDV2T, U1, LDU1, 410 $ U2, LDU2, WORK, LWORK, IWORK, INFO ) 411 RETURN 412 END IF 413* 414* Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if 415* convenient 416* 417 IF( INFO .EQ. 0 .AND. M-Q .LT. Q ) THEN 418 IF( DEFAULTSIGNS ) THEN 419 SIGNST = 'O' 420 ELSE 421 SIGNST = 'D' 422 END IF 423 CALL DORCSD( JOBU2, JOBU1, JOBV2T, JOBV1T, TRANS, SIGNST, M, 424 $ M-P, M-Q, X22, LDX22, X21, LDX21, X12, LDX12, X11, 425 $ LDX11, THETA, U2, LDU2, U1, LDU1, V2T, LDV2T, V1T, 426 $ LDV1T, WORK, LWORK, IWORK, INFO ) 427 RETURN 428 END IF 429* 430* Compute workspace 431* 432 IF( INFO .EQ. 0 ) THEN 433* 434 IPHI = 2 435 ITAUP1 = IPHI + MAX( 1, Q - 1 ) 436 ITAUP2 = ITAUP1 + MAX( 1, P ) 437 ITAUQ1 = ITAUP2 + MAX( 1, M - P ) 438 ITAUQ2 = ITAUQ1 + MAX( 1, Q ) 439 IORGQR = ITAUQ2 + MAX( 1, M - Q ) 440 CALL DORGQR( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 441 $ CHILDINFO ) 442 LORGQRWORKOPT = INT( WORK(1) ) 443 LORGQRWORKMIN = MAX( 1, M - Q ) 444 IORGLQ = ITAUQ2 + MAX( 1, M - Q ) 445 CALL DORGLQ( M-Q, M-Q, M-Q, U1, MAX(1,M-Q), U1, WORK, -1, 446 $ CHILDINFO ) 447 LORGLQWORKOPT = INT( WORK(1) ) 448 LORGLQWORKMIN = MAX( 1, M - Q ) 449 IORBDB = ITAUQ2 + MAX( 1, M - Q ) 450 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, 451 $ X21, LDX21, X22, LDX22, THETA, V1T, U1, U2, V1T, 452 $ V2T, WORK, -1, CHILDINFO ) 453 LORBDBWORKOPT = INT( WORK(1) ) 454 LORBDBWORKMIN = LORBDBWORKOPT 455 IB11D = ITAUQ2 + MAX( 1, M - Q ) 456 IB11E = IB11D + MAX( 1, Q ) 457 IB12D = IB11E + MAX( 1, Q - 1 ) 458 IB12E = IB12D + MAX( 1, Q ) 459 IB21D = IB12E + MAX( 1, Q - 1 ) 460 IB21E = IB21D + MAX( 1, Q ) 461 IB22D = IB21E + MAX( 1, Q - 1 ) 462 IB22E = IB22D + MAX( 1, Q ) 463 IBBCSD = IB22E + MAX( 1, Q - 1 ) 464 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, 465 $ THETA, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 466 $ LDV2T, U1, U1, U1, U1, U1, U1, U1, U1, WORK, -1, 467 $ CHILDINFO ) 468 LBBCSDWORKOPT = INT( WORK(1) ) 469 LBBCSDWORKMIN = LBBCSDWORKOPT 470 LWORKOPT = MAX( IORGQR + LORGQRWORKOPT, IORGLQ + LORGLQWORKOPT, 471 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKOPT ) - 1 472 LWORKMIN = MAX( IORGQR + LORGQRWORKMIN, IORGLQ + LORGLQWORKMIN, 473 $ IORBDB + LORBDBWORKOPT, IBBCSD + LBBCSDWORKMIN ) - 1 474 WORK(1) = MAX(LWORKOPT,LWORKMIN) 475* 476 IF( LWORK .LT. LWORKMIN .AND. .NOT. LQUERY ) THEN 477 INFO = -22 478 ELSE 479 LORGQRWORK = LWORK - IORGQR + 1 480 LORGLQWORK = LWORK - IORGLQ + 1 481 LORBDBWORK = LWORK - IORBDB + 1 482 LBBCSDWORK = LWORK - IBBCSD + 1 483 END IF 484 END IF 485* 486* Abort if any illegal arguments 487* 488 IF( INFO .NE. 0 ) THEN 489 CALL XERBLA( 'DORCSD', -INFO ) 490 RETURN 491 ELSE IF( LQUERY ) THEN 492 RETURN 493 END IF 494* 495* Transform to bidiagonal block form 496* 497 CALL DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, 498 $ LDX21, X22, LDX22, THETA, WORK(IPHI), WORK(ITAUP1), 499 $ WORK(ITAUP2), WORK(ITAUQ1), WORK(ITAUQ2), 500 $ WORK(IORBDB), LORBDBWORK, CHILDINFO ) 501* 502* Accumulate Householder reflectors 503* 504 IF( COLMAJOR ) THEN 505 IF( WANTU1 .AND. P .GT. 0 ) THEN 506 CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 ) 507 CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR), 508 $ LORGQRWORK, INFO) 509 END IF 510 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 511 CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 ) 512 CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 513 $ WORK(IORGQR), LORGQRWORK, INFO ) 514 END IF 515 IF( WANTV1T .AND. Q .GT. 0 ) THEN 516 CALL DLACPY( 'U', Q-1, Q-1, X11(1,2), LDX11, V1T(2,2), 517 $ LDV1T ) 518 V1T(1, 1) = ONE 519 DO J = 2, Q 520 V1T(1,J) = ZERO 521 V1T(J,1) = ZERO 522 END DO 523 CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 524 $ WORK(IORGLQ), LORGLQWORK, INFO ) 525 END IF 526 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 527 CALL DLACPY( 'U', P, M-Q, X12, LDX12, V2T, LDV2T ) 528 IF (M-P .GT. Q) Then 529 CALL DLACPY( 'U', M-P-Q, M-P-Q, X22(Q+1,P+1), LDX22, 530 $ V2T(P+1,P+1), LDV2T ) 531 END IF 532 IF (M .GT. Q) THEN 533 CALL DORGLQ( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 534 $ WORK(IORGLQ), LORGLQWORK, INFO ) 535 END IF 536 END IF 537 ELSE 538 IF( WANTU1 .AND. P .GT. 0 ) THEN 539 CALL DLACPY( 'U', Q, P, X11, LDX11, U1, LDU1 ) 540 CALL DORGLQ( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGLQ), 541 $ LORGLQWORK, INFO) 542 END IF 543 IF( WANTU2 .AND. M-P .GT. 0 ) THEN 544 CALL DLACPY( 'U', Q, M-P, X21, LDX21, U2, LDU2 ) 545 CALL DORGLQ( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2), 546 $ WORK(IORGLQ), LORGLQWORK, INFO ) 547 END IF 548 IF( WANTV1T .AND. Q .GT. 0 ) THEN 549 CALL DLACPY( 'L', Q-1, Q-1, X11(2,1), LDX11, V1T(2,2), 550 $ LDV1T ) 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 DORGQR( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1), 557 $ WORK(IORGQR), LORGQRWORK, INFO ) 558 END IF 559 IF( WANTV2T .AND. M-Q .GT. 0 ) THEN 560 CALL DLACPY( 'L', M-Q, P, X12, LDX12, V2T, LDV2T ) 561 CALL DLACPY( 'L', M-P-Q, M-P-Q, X22(P+1,Q+1), LDX22, 562 $ V2T(P+1,P+1), LDV2T ) 563 CALL DORGQR( M-Q, M-Q, M-Q, V2T, LDV2T, WORK(ITAUQ2), 564 $ WORK(IORGQR), LORGQRWORK, INFO ) 565 END IF 566 END IF 567* 568* Compute the CSD of the matrix in bidiagonal-block form 569* 570 CALL DBBCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, 571 $ WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, 572 $ LDV2T, WORK(IB11D), WORK(IB11E), WORK(IB12D), 573 $ WORK(IB12E), WORK(IB21D), WORK(IB21E), WORK(IB22D), 574 $ WORK(IB22E), WORK(IBBCSD), LBBCSDWORK, INFO ) 575* 576* Permute rows and columns to place identity submatrices in top- 577* left corner of (1,1)-block and/or bottom-right corner of (1,2)- 578* block and/or bottom-right corner of (2,1)-block and/or top-left 579* corner of (2,2)-block 580* 581 IF( Q .GT. 0 .AND. WANTU2 ) THEN 582 DO I = 1, Q 583 IWORK(I) = M - P - Q + I 584 END DO 585 DO I = Q + 1, M - P 586 IWORK(I) = I - Q 587 END DO 588 IF( COLMAJOR ) THEN 589 CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 590 ELSE 591 CALL DLAPMR( .FALSE., M-P, M-P, U2, LDU2, IWORK ) 592 END IF 593 END IF 594 IF( M .GT. 0 .AND. WANTV2T ) THEN 595 DO I = 1, P 596 IWORK(I) = M - P - Q + I 597 END DO 598 DO I = P + 1, M - Q 599 IWORK(I) = I - P 600 END DO 601 IF( .NOT. COLMAJOR ) THEN 602 CALL DLAPMT( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 603 ELSE 604 CALL DLAPMR( .FALSE., M-Q, M-Q, V2T, LDV2T, IWORK ) 605 END IF 606 END IF 607* 608 RETURN 609* 610* End DORCSD 611* 612 END 613 614