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