LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zbbcsd()

subroutine zbbcsd ( character jobu1,
character jobu2,
character jobv1t,
character jobv2t,
character trans,
integer m,
integer p,
integer q,
double precision, dimension( * ) theta,
double precision, dimension( * ) phi,
complex*16, dimension( ldu1, * ) u1,
integer ldu1,
complex*16, dimension( ldu2, * ) u2,
integer ldu2,
complex*16, dimension( ldv1t, * ) v1t,
integer ldv1t,
complex*16, dimension( ldv2t, * ) v2t,
integer ldv2t,
double precision, dimension( * ) b11d,
double precision, dimension( * ) b11e,
double precision, dimension( * ) b12d,
double precision, dimension( * ) b12e,
double precision, dimension( * ) b21d,
double precision, dimension( * ) b21e,
double precision, dimension( * ) b22d,
double precision, dimension( * ) b22e,
double precision, dimension( * ) rwork,
integer lrwork,
integer info )

ZBBCSD

Download ZBBCSD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZBBCSD computes the CS decomposition of a unitary matrix in
!> bidiagonal-block form,
!>
!>
!>     [ B11 | B12 0  0 ]
!>     [  0  |  0 -I  0 ]
!> X = [----------------]
!>     [ B21 | B22 0  0 ]
!>     [  0  |  0  0  I ]
!>
!>                               [  C | -S  0  0 ]
!>                   [ U1 |    ] [  0 |  0 -I  0 ] [ V1 |    ]**H
!>                 = [---------] [---------------] [---------]   .
!>                   [    | U2 ] [  S |  C  0  0 ] [    | V2 ]
!>                               [  0 |  0  0  I ]
!>
!> X is M-by-M, its top-left block is P-by-Q, and Q must be no larger
!> than P, M-P, or M-Q. (If Q is not the smallest index, then X must be
!> transposed and/or permuted. This can be done in constant time using
!> the TRANS and SIGNS options. See ZUNCSD for details.)
!>
!> The bidiagonal matrices B11, B12, B21, and B22 are represented
!> implicitly by angles THETA(1:Q) and PHI(1:Q-1).
!>
!> The unitary matrices U1, U2, V1T, and V2T are input/output.
!> The input matrices are pre- or post-multiplied by the appropriate
!> singular vector matrices.
!>
Parameters
[in]JOBU1
!>          JOBU1 is CHARACTER
!>          = 'Y':      U1 is updated;
!>          otherwise:  U1 is not updated.
!>
[in]JOBU2
!>          JOBU2 is CHARACTER
!>          = 'Y':      U2 is updated;
!>          otherwise:  U2 is not updated.
!>
[in]JOBV1T
!>          JOBV1T is CHARACTER
!>          = 'Y':      V1T is updated;
!>          otherwise:  V1T is not updated.
!>
[in]JOBV2T
!>          JOBV2T is CHARACTER
!>          = 'Y':      V2T is updated;
!>          otherwise:  V2T is not updated.
!>
[in]TRANS
!>          TRANS is CHARACTER
!>          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
!>                      order;
!>          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
!>                      major order.
!>
[in]M
!>          M is INTEGER
!>          The number of rows and columns in X, the unitary matrix in
!>          bidiagonal-block form.
!>
[in]P
!>          P is INTEGER
!>          The number of rows in the top-left block of X. 0 <= P <= M.
!>
[in]Q
!>          Q is INTEGER
!>          The number of columns in the top-left block of X.
!>          0 <= Q <= MIN(P,M-P,M-Q).
!>
[in,out]THETA
!>          THETA is DOUBLE PRECISION array, dimension (Q)
!>          On entry, the angles THETA(1),...,THETA(Q) that, along with
!>          PHI(1), ...,PHI(Q-1), define the matrix in bidiagonal-block
!>          form. On exit, the angles whose cosines and sines define the
!>          diagonal blocks in the CS decomposition.
!>
[in,out]PHI
!>          PHI is DOUBLE PRECISION array, dimension (Q-1)
!>          The angles PHI(1),...,PHI(Q-1) that, along with THETA(1),...,
!>          THETA(Q), define the matrix in bidiagonal-block form.
!>
[in,out]U1
!>          U1 is COMPLEX*16 array, dimension (LDU1,P)
!>          On entry, a P-by-P matrix. On exit, U1 is postmultiplied
!>          by the left singular vector matrix common to [ B11 ; 0 ] and
!>          [ B12 0 0 ; 0 -I 0 0 ].
!>
[in]LDU1
!>          LDU1 is INTEGER
!>          The leading dimension of the array U1, LDU1 >= MAX(1,P).
!>
[in,out]U2
!>          U2 is COMPLEX*16 array, dimension (LDU2,M-P)
!>          On entry, an (M-P)-by-(M-P) matrix. On exit, U2 is
!>          postmultiplied by the left singular vector matrix common to
!>          [ B21 ; 0 ] and [ B22 0 0 ; 0 0 I ].
!>
[in]LDU2
!>          LDU2 is INTEGER
!>          The leading dimension of the array U2, LDU2 >= MAX(1,M-P).
!>
[in,out]V1T
!>          V1T is COMPLEX*16 array, dimension (LDV1T,Q)
!>          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
!>          by the conjugate transpose of the right singular vector
!>          matrix common to [ B11 ; 0 ] and [ B21 ; 0 ].
!>
[in]LDV1T
!>          LDV1T is INTEGER
!>          The leading dimension of the array V1T, LDV1T >= MAX(1,Q).
!>
[in,out]V2T
!>          V2T is COMPLEX*16 array, dimension (LDV2T,M-Q)
!>          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
!>          premultiplied by the conjugate transpose of the right
!>          singular vector matrix common to [ B12 0 0 ; 0 -I 0 ] and
!>          [ B22 0 0 ; 0 0 I ].
!>
[in]LDV2T
!>          LDV2T is INTEGER
!>          The leading dimension of the array V2T, LDV2T >= MAX(1,M-Q).
!>
[out]B11D
!>          B11D is DOUBLE PRECISION array, dimension (Q)
!>          When ZBBCSD converges, B11D contains the cosines of THETA(1),
!>          ..., THETA(Q). If ZBBCSD fails to converge, then B11D
!>          contains the diagonal of the partially reduced top-left
!>          block.
!>
[out]B11E
!>          B11E is DOUBLE PRECISION array, dimension (Q-1)
!>          When ZBBCSD converges, B11E contains zeros. If ZBBCSD fails
!>          to converge, then B11E contains the superdiagonal of the
!>          partially reduced top-left block.
!>
[out]B12D
!>          B12D is DOUBLE PRECISION array, dimension (Q)
!>          When ZBBCSD converges, B12D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
!>          B12D contains the diagonal of the partially reduced top-right
!>          block.
!>
[out]B12E
!>          B12E is DOUBLE PRECISION array, dimension (Q-1)
!>          When ZBBCSD converges, B12E contains zeros. If ZBBCSD fails
!>          to converge, then B12E contains the subdiagonal of the
!>          partially reduced top-right block.
!>
[out]B21D
!>          B21D is DOUBLE PRECISION array, dimension (Q)
!>          When ZBBCSD converges, B21D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
!>          B21D contains the diagonal of the partially reduced bottom-left
!>          block.
!>
[out]B21E
!>          B21E is DOUBLE PRECISION array, dimension (Q-1)
!>          When ZBBCSD converges, B21E contains zeros. If ZBBCSD fails
!>          to converge, then B21E contains the subdiagonal of the
!>          partially reduced bottom-left block.
!>
[out]B22D
!>          B22D is DOUBLE PRECISION array, dimension (Q)
!>          When ZBBCSD converges, B22D contains the negative sines of
!>          THETA(1), ..., THETA(Q). If ZBBCSD fails to converge, then
!>          B22D contains the diagonal of the partially reduced bottom-right
!>          block.
!>
[out]B22E
!>          B22E is DOUBLE PRECISION array, dimension (Q-1)
!>          When ZBBCSD converges, B22E contains zeros. If ZBBCSD fails
!>          to converge, then B22E contains the subdiagonal of the
!>          partially reduced bottom-right block.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK. LRWORK >= MAX(1,8*Q).
!>
!>          If LRWORK = -1, then a workspace query is assumed; the
!>          routine only calculates the optimal size of the RWORK array,
!>          returns this value as the first entry of the work array, and
!>          no error message related to LRWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  if ZBBCSD did not converge, INFO specifies the number
!>                of nonzero entries in PHI, and B11D, B11E, etc.,
!>                contain the partially reduced matrix.
!>
Internal Parameters:
!>  TOLMUL  DOUBLE PRECISION, default = MAX(10,MIN(100,EPS**(-1/8)))
!>          TOLMUL controls the convergence criterion of the QR loop.
!>          Angles THETA(i), PHI(i) are rounded to 0 or PI/2 when they
!>          are within TOLMUL*EPS of either bound.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 326 of file zbbcsd.f.

331*
332* -- LAPACK computational routine --
333* -- LAPACK is a software package provided by Univ. of Tennessee, --
334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335*
336* .. Scalar Arguments ..
337 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
338 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
339* ..
340* .. Array Arguments ..
341 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
342 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
343 $ PHI( * ), THETA( * ), RWORK( * )
344 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
345 $ V2T( LDV2T, * )
346* ..
347*
348* ===================================================================
349*
350* .. Parameters ..
351 INTEGER MAXITR
352 parameter( maxitr = 6 )
353 DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, TEN, ZERO
354 parameter( hundred = 100.0d0, meighth = -0.125d0,
355 $ one = 1.0d0, ten = 10.0d0, zero = 0.0d0 )
356 COMPLEX*16 NEGONECOMPLEX
357 parameter( negonecomplex = (-1.0d0,0.0d0) )
358 DOUBLE PRECISION PIOVER2
359 parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
360* ..
361* .. Local Scalars ..
362 LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
363 $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
364 $ WANTV2T
365 INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
366 $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
367 $ LRWORKMIN, LRWORKOPT, MAXIT, MINI
368 DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
369 $ EPS, MU, NU, R, SIGMA11, SIGMA21,
370 $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
371 $ UNFL, X1, X2, Y1, Y2
372*
373 EXTERNAL dlartgp, dlartgs, dlas2, xerbla, zlasr,
374 $ zscal,
375 $ zswap
376* ..
377* .. External Functions ..
378 DOUBLE PRECISION DLAMCH
379 LOGICAL LSAME
380 EXTERNAL lsame, dlamch
381* ..
382* .. Intrinsic Functions ..
383 INTRINSIC abs, atan2, cos, max, min, sin, sqrt
384* ..
385* .. Executable Statements ..
386*
387* Test input arguments
388*
389 info = 0
390 lquery = lrwork .EQ. -1
391 wantu1 = lsame( jobu1, 'Y' )
392 wantu2 = lsame( jobu2, 'Y' )
393 wantv1t = lsame( jobv1t, 'Y' )
394 wantv2t = lsame( jobv2t, 'Y' )
395 colmajor = .NOT. lsame( trans, 'T' )
396*
397 IF( m .LT. 0 ) THEN
398 info = -6
399 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
400 info = -7
401 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
402 info = -8
403 ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q ) THEN
404 info = -8
405 ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
406 info = -12
407 ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
408 info = -14
409 ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
410 info = -16
411 ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
412 info = -18
413 END IF
414*
415* Quick return if Q = 0
416*
417 IF( info .EQ. 0 .AND. q .EQ. 0 ) THEN
418 lrworkmin = 1
419 rwork(1) = lrworkmin
420 RETURN
421 END IF
422*
423* Compute workspace
424*
425 IF( info .EQ. 0 ) THEN
426 iu1cs = 1
427 iu1sn = iu1cs + q
428 iu2cs = iu1sn + q
429 iu2sn = iu2cs + q
430 iv1tcs = iu2sn + q
431 iv1tsn = iv1tcs + q
432 iv2tcs = iv1tsn + q
433 iv2tsn = iv2tcs + q
434 lrworkopt = iv2tsn + q - 1
435 lrworkmin = lrworkopt
436 rwork(1) = lrworkopt
437 IF( lrwork .LT. lrworkmin .AND. .NOT. lquery ) THEN
438 info = -28
439 END IF
440 END IF
441*
442 IF( info .NE. 0 ) THEN
443 CALL xerbla( 'ZBBCSD', -info )
444 RETURN
445 ELSE IF( lquery ) THEN
446 RETURN
447 END IF
448*
449* Get machine constants
450*
451 eps = dlamch( 'Epsilon' )
452 unfl = dlamch( 'Safe minimum' )
453 tolmul = max( ten, min( hundred, eps**meighth ) )
454 tol = tolmul*eps
455 thresh = max( tol, maxitr*q*q*unfl )
456*
457* Test for negligible sines or cosines
458*
459 DO i = 1, q
460 IF( theta(i) .LT. thresh ) THEN
461 theta(i) = zero
462 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
463 theta(i) = piover2
464 END IF
465 END DO
466 DO i = 1, q-1
467 IF( phi(i) .LT. thresh ) THEN
468 phi(i) = zero
469 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
470 phi(i) = piover2
471 END IF
472 END DO
473*
474* Initial deflation
475*
476 imax = q
477 DO WHILE( imax .GT. 1 )
478 IF( phi(imax-1) .NE. zero ) THEN
479 EXIT
480 END IF
481 imax = imax - 1
482 END DO
483 imin = imax - 1
484 IF ( imin .GT. 1 ) THEN
485 DO WHILE( phi(imin-1) .NE. zero )
486 imin = imin - 1
487 IF ( imin .LE. 1 ) EXIT
488 END DO
489 END IF
490*
491* Initialize iteration counter
492*
493 maxit = maxitr*q*q
494 iter = 0
495*
496* Begin main iteration loop
497*
498 DO WHILE( imax .GT. 1 )
499*
500* Compute the matrix entries
501*
502 b11d(imin) = cos( theta(imin) )
503 b21d(imin) = -sin( theta(imin) )
504 DO i = imin, imax - 1
505 b11e(i) = -sin( theta(i) ) * sin( phi(i) )
506 b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
507 b12d(i) = sin( theta(i) ) * cos( phi(i) )
508 b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
509 b21e(i) = -cos( theta(i) ) * sin( phi(i) )
510 b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
511 b22d(i) = cos( theta(i) ) * cos( phi(i) )
512 b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
513 END DO
514 b12d(imax) = sin( theta(imax) )
515 b22d(imax) = cos( theta(imax) )
516*
517* Abort if not converging; otherwise, increment ITER
518*
519 IF( iter .GT. maxit ) THEN
520 info = 0
521 DO i = 1, q
522 IF( phi(i) .NE. zero )
523 $ info = info + 1
524 END DO
525 RETURN
526 END IF
527*
528 iter = iter + imax - imin
529*
530* Compute shifts
531*
532 thetamax = theta(imin)
533 thetamin = theta(imin)
534 DO i = imin+1, imax
535 IF( theta(i) > thetamax )
536 $ thetamax = theta(i)
537 IF( theta(i) < thetamin )
538 $ thetamin = theta(i)
539 END DO
540*
541 IF( thetamax .GT. piover2 - thresh ) THEN
542*
543* Zero on diagonals of B11 and B22; induce deflation with a
544* zero shift
545*
546 mu = zero
547 nu = one
548*
549 ELSE IF( thetamin .LT. thresh ) THEN
550*
551* Zero on diagonals of B12 and B22; induce deflation with a
552* zero shift
553*
554 mu = one
555 nu = zero
556*
557 ELSE
558*
559* Compute shifts for B11 and B21 and use the lesser
560*
561 CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax),
562 $ sigma11,
563 $ dummy )
564 CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax),
565 $ sigma21,
566 $ dummy )
567*
568 IF( sigma11 .LE. sigma21 ) THEN
569 mu = sigma11
570 nu = sqrt( one - mu**2 )
571 IF( mu .LT. thresh ) THEN
572 mu = zero
573 nu = one
574 END IF
575 ELSE
576 nu = sigma21
577 mu = sqrt( 1.0 - nu**2 )
578 IF( nu .LT. thresh ) THEN
579 mu = one
580 nu = zero
581 END IF
582 END IF
583 END IF
584*
585* Rotate to produce bulges in B11 and B21
586*
587 IF( mu .LE. nu ) THEN
588 CALL dlartgs( b11d(imin), b11e(imin), mu,
589 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
590 ELSE
591 CALL dlartgs( b21d(imin), b21e(imin), nu,
592 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
593 END IF
594*
595 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
596 $ rwork(iv1tsn+imin-1)*b11e(imin)
597 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
598 $ rwork(iv1tsn+imin-1)*b11d(imin)
599 b11d(imin) = temp
600 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
601 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
602 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
603 $ rwork(iv1tsn+imin-1)*b21e(imin)
604 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
605 $ rwork(iv1tsn+imin-1)*b21d(imin)
606 b21d(imin) = temp
607 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
608 b21d(imin+1) = rwork(iv1tcs+imin-1)*b21d(imin+1)
609*
610* Compute THETA(IMIN)
611*
612 theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
613 $ sqrt( b11d(imin)**2+b11bulge**2 ) )
614*
615* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
616*
617 IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 ) THEN
618 CALL dlartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
619 $ rwork(iu1cs+imin-1), r )
620 ELSE IF( mu .LE. nu ) THEN
621 CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
622 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
623 ELSE
624 CALL dlartgs( b12d( imin ), b12e( imin ), nu,
625 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
626 END IF
627 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
628 CALL dlartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
629 $ rwork(iu2cs+imin-1), r )
630 ELSE IF( nu .LT. mu ) THEN
631 CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
632 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
633 ELSE
634 CALL dlartgs( b22d(imin), b22e(imin), mu,
635 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
636 END IF
637 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
638 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
639*
640 temp = rwork(iu1cs+imin-1)*b11e(imin) +
641 $ rwork(iu1sn+imin-1)*b11d(imin+1)
642 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
643 $ rwork(iu1sn+imin-1)*b11e(imin)
644 b11e(imin) = temp
645 IF( imax .GT. imin+1 ) THEN
646 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
647 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
648 END IF
649 temp = rwork(iu1cs+imin-1)*b12d(imin) +
650 $ rwork(iu1sn+imin-1)*b12e(imin)
651 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
652 $ rwork(iu1sn+imin-1)*b12d(imin)
653 b12d(imin) = temp
654 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
655 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
656 temp = rwork(iu2cs+imin-1)*b21e(imin) +
657 $ rwork(iu2sn+imin-1)*b21d(imin+1)
658 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
659 $ rwork(iu2sn+imin-1)*b21e(imin)
660 b21e(imin) = temp
661 IF( imax .GT. imin+1 ) THEN
662 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
663 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
664 END IF
665 temp = rwork(iu2cs+imin-1)*b22d(imin) +
666 $ rwork(iu2sn+imin-1)*b22e(imin)
667 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
668 $ rwork(iu2sn+imin-1)*b22d(imin)
669 b22d(imin) = temp
670 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
671 b22d(imin+1) = rwork(iu2cs+imin-1)*b22d(imin+1)
672*
673* Inner loop: chase bulges from B11(IMIN,IMIN+2),
674* B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
675* bottom-right
676*
677 DO i = imin+1, imax-1
678*
679* Compute PHI(I-1)
680*
681 x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
682 x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
683 y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
684 y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
685*
686 phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
687*
688* Determine if there are bulges to chase or if a new direct
689* summand has been reached
690*
691 restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
692 restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
693 restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
694 restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
695*
696* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
697* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
698* chasing by applying the original shift again.
699*
700 IF( .NOT. restart11 .AND. .NOT. restart21 ) THEN
701 CALL dlartgp( x2, x1, rwork(iv1tsn+i-1),
702 $ rwork(iv1tcs+i-1), r )
703 ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
704 CALL dlartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
705 $ rwork(iv1tcs+i-1), r )
706 ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
707 CALL dlartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
708 $ rwork(iv1tcs+i-1), r )
709 ELSE IF( mu .LE. nu ) THEN
710 CALL dlartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
711 $ rwork(iv1tsn+i-1) )
712 ELSE
713 CALL dlartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
714 $ rwork(iv1tsn+i-1) )
715 END IF
716 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
717 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
718 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
719 CALL dlartgp( y2, y1, rwork(iv2tsn+i-1-1),
720 $ rwork(iv2tcs+i-1-1), r )
721 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
722 CALL dlartgp( b12bulge, b12d(i-1),
723 $ rwork(iv2tsn+i-1-1),
724 $ rwork(iv2tcs+i-1-1), r )
725 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
726 CALL dlartgp( b22bulge, b22d(i-1),
727 $ rwork(iv2tsn+i-1-1),
728 $ rwork(iv2tcs+i-1-1), r )
729 ELSE IF( nu .LT. mu ) THEN
730 CALL dlartgs( b12e(i-1), b12d(i), nu,
731 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
732 ELSE
733 CALL dlartgs( b22e(i-1), b22d(i), mu,
734 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
735 END IF
736*
737 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
738 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
739 $ rwork(iv1tsn+i-1)*b11d(i)
740 b11d(i) = temp
741 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
742 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
743 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
744 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
745 $ rwork(iv1tsn+i-1)*b21d(i)
746 b21d(i) = temp
747 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
748 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
749 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
750 $ rwork(iv2tsn+i-1-1)*b12d(i)
751 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
752 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
753 b12e(i-1) = temp
754 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
755 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
756 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
757 $ rwork(iv2tsn+i-1-1)*b22d(i)
758 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
759 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
760 b22e(i-1) = temp
761 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
762 b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
763*
764* Compute THETA(I)
765*
766 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
767 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
768 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
769 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
770*
771 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
772*
773* Determine if there are bulges to chase or if a new direct
774* summand has been reached
775*
776 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
777 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
778 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
779 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
780*
781* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
782* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
783* chasing by applying the original shift again.
784*
785 IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
786 CALL dlartgp( x2, x1, rwork(iu1sn+i-1),
787 $ rwork(iu1cs+i-1),
788 $ r )
789 ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
790 CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
791 $ rwork(iu1cs+i-1), r )
792 ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
793 CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
794 $ rwork(iu1cs+i-1), r )
795 ELSE IF( mu .LE. nu ) THEN
796 CALL dlartgs( b11e(i), b11d(i+1), mu,
797 $ rwork(iu1cs+i-1),
798 $ rwork(iu1sn+i-1) )
799 ELSE
800 CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
801 $ rwork(iu1sn+i-1) )
802 END IF
803 IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
804 CALL dlartgp( y2, y1, rwork(iu2sn+i-1),
805 $ rwork(iu2cs+i-1),
806 $ r )
807 ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
808 CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
809 $ rwork(iu2cs+i-1), r )
810 ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
811 CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
812 $ rwork(iu2cs+i-1), r )
813 ELSE IF( nu .LT. mu ) THEN
814 CALL dlartgs( b21e(i), b21e(i+1), nu,
815 $ rwork(iu2cs+i-1),
816 $ rwork(iu2sn+i-1) )
817 ELSE
818 CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
819 $ rwork(iu2sn+i-1) )
820 END IF
821 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
822 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
823*
824 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
825 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
826 $ rwork(iu1sn+i-1)*b11e(i)
827 b11e(i) = temp
828 IF( i .LT. imax - 1 ) THEN
829 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
830 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
831 END IF
832 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
833 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
834 $ rwork(iu2sn+i-1)*b21e(i)
835 b21e(i) = temp
836 IF( i .LT. imax - 1 ) THEN
837 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
838 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
839 END IF
840 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
841 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
842 $ rwork(iu1sn+i-1)*b12d(i)
843 b12d(i) = temp
844 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
845 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
846 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
847 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
848 $ rwork(iu2sn+i-1)*b22d(i)
849 b22d(i) = temp
850 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
851 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
852*
853 END DO
854*
855* Compute PHI(IMAX-1)
856*
857 x1 = sin(theta(imax-1))*b11e(imax-1) +
858 $ cos(theta(imax-1))*b21e(imax-1)
859 y1 = sin(theta(imax-1))*b12d(imax-1) +
860 $ cos(theta(imax-1))*b22d(imax-1)
861 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
862*
863 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
864*
865* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
866*
867 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
868 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
869*
870 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
871 CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
872 $ rwork(iv2tcs+imax-1-1), r )
873 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
874 CALL dlartgp( b12bulge, b12d(imax-1),
875 $ rwork(iv2tsn+imax-1-1),
876 $ rwork(iv2tcs+imax-1-1), r )
877 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
878 CALL dlartgp( b22bulge, b22d(imax-1),
879 $ rwork(iv2tsn+imax-1-1),
880 $ rwork(iv2tcs+imax-1-1), r )
881 ELSE IF( nu .LT. mu ) THEN
882 CALL dlartgs( b12e(imax-1), b12d(imax), nu,
883 $ rwork(iv2tcs+imax-1-1),
884 $ rwork(iv2tsn+imax-1-1) )
885 ELSE
886 CALL dlartgs( b22e(imax-1), b22d(imax), mu,
887 $ rwork(iv2tcs+imax-1-1),
888 $ rwork(iv2tsn+imax-1-1) )
889 END IF
890*
891 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
892 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
893 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
894 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
895 b12e(imax-1) = temp
896 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
897 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
898 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
899 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
900 b22e(imax-1) = temp
901*
902* Update singular vectors
903*
904 IF( wantu1 ) THEN
905 IF( colmajor ) THEN
906 CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
907 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
908 $ u1(1,imin), ldu1 )
909 ELSE
910 CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
911 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
912 $ u1(imin,1), ldu1 )
913 END IF
914 END IF
915 IF( wantu2 ) THEN
916 IF( colmajor ) THEN
917 CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
918 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
919 $ u2(1,imin), ldu2 )
920 ELSE
921 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
922 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
923 $ u2(imin,1), ldu2 )
924 END IF
925 END IF
926 IF( wantv1t ) THEN
927 IF( colmajor ) THEN
928 CALL zlasr( 'L', 'V', 'F', imax-imin+1, q,
929 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
930 $ v1t(imin,1), ldv1t )
931 ELSE
932 CALL zlasr( 'R', 'V', 'F', q, imax-imin+1,
933 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
934 $ v1t(1,imin), ldv1t )
935 END IF
936 END IF
937 IF( wantv2t ) THEN
938 IF( colmajor ) THEN
939 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-q,
940 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
941 $ v2t(imin,1), ldv2t )
942 ELSE
943 CALL zlasr( 'R', 'V', 'F', m-q, imax-imin+1,
944 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
945 $ v2t(1,imin), ldv2t )
946 END IF
947 END IF
948*
949* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
950*
951 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
952 b11d(imax) = -b11d(imax)
953 b21d(imax) = -b21d(imax)
954 IF( wantv1t ) THEN
955 IF( colmajor ) THEN
956 CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
957 ELSE
958 CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
959 END IF
960 END IF
961 END IF
962*
963* Compute THETA(IMAX)
964*
965 x1 = cos(phi(imax-1))*b11d(imax) +
966 $ sin(phi(imax-1))*b12e(imax-1)
967 y1 = cos(phi(imax-1))*b21d(imax) +
968 $ sin(phi(imax-1))*b22e(imax-1)
969*
970 theta(imax) = atan2( abs(y1), abs(x1) )
971*
972* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
973* and B22(IMAX,IMAX-1)
974*
975 IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
976 b12d(imax) = -b12d(imax)
977 IF( wantu1 ) THEN
978 IF( colmajor ) THEN
979 CALL zscal( p, negonecomplex, u1(1,imax), 1 )
980 ELSE
981 CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
982 END IF
983 END IF
984 END IF
985 IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
986 b22d(imax) = -b22d(imax)
987 IF( wantu2 ) THEN
988 IF( colmajor ) THEN
989 CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
990 ELSE
991 CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
992 END IF
993 END IF
994 END IF
995*
996* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
997*
998 IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
999 IF( wantv2t ) THEN
1000 IF( colmajor ) THEN
1001 CALL zscal( m-q, negonecomplex, v2t(imax,1),
1002 $ ldv2t )
1003 ELSE
1004 CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
1005 END IF
1006 END IF
1007 END IF
1008*
1009* Test for negligible sines or cosines
1010*
1011 DO i = imin, imax
1012 IF( theta(i) .LT. thresh ) THEN
1013 theta(i) = zero
1014 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1015 theta(i) = piover2
1016 END IF
1017 END DO
1018 DO i = imin, imax-1
1019 IF( phi(i) .LT. thresh ) THEN
1020 phi(i) = zero
1021 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1022 phi(i) = piover2
1023 END IF
1024 END DO
1025*
1026* Deflate
1027*
1028 IF (imax .GT. 1) THEN
1029 DO WHILE( phi(imax-1) .EQ. zero )
1030 imax = imax - 1
1031 IF (imax .LE. 1) EXIT
1032 END DO
1033 END IF
1034 IF( imin .GT. imax - 1 )
1035 $ imin = imax - 1
1036 IF (imin .GT. 1) THEN
1037 DO WHILE (phi(imin-1) .NE. zero)
1038 imin = imin - 1
1039 IF (imin .LE. 1) EXIT
1040 END DO
1041 END IF
1042*
1043* Repeat main iteration loop
1044*
1045 END DO
1046*
1047* Postprocessing: order THETA from least to greatest
1048*
1049 DO i = 1, q
1050*
1051 mini = i
1052 thetamin = theta(i)
1053 DO j = i+1, q
1054 IF( theta(j) .LT. thetamin ) THEN
1055 mini = j
1056 thetamin = theta(j)
1057 END IF
1058 END DO
1059*
1060 IF( mini .NE. i ) THEN
1061 theta(mini) = theta(i)
1062 theta(i) = thetamin
1063 IF( colmajor ) THEN
1064 IF( wantu1 )
1065 $ CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1066 IF( wantu2 )
1067 $ CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1068 IF( wantv1t )
1069 $ CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1),
1070 $ ldv1t )
1071 IF( wantv2t )
1072 $ CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1073 $ ldv2t )
1074 ELSE
1075 IF( wantu1 )
1076 $ CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1077 IF( wantu2 )
1078 $ CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1079 IF( wantv1t )
1080 $ CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1081 IF( wantv2t )
1082 $ CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1083 END IF
1084 END IF
1085*
1086 END DO
1087*
1088 RETURN
1089*
1090* End of ZBBCSD
1091*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlartgp
subroutine dlartgp(f, g, cs, sn, r)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition dlartgp.f:93
dlartgs
subroutine dlartgs(x, y, sigma, cs, sn)
DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
Definition dlartgs.f:88
dlas2
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:103
zlasr
subroutine zlasr(side, pivot, direct, m, n, c, s, a, lda)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition zlasr.f:198
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zscal
subroutine zscal(n, za, zx, incx)
ZSCAL
Definition zscal.f:78
zswap
subroutine zswap(n, zx, incx, zy, incy)
ZSWAP
Definition zswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: