LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, dimenison (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,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[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.
Date
June 2016

Definition at line 334 of file zbbcsd.f.

334 *
335 * -- LAPACK computational routine (version 3.6.1) --
336 * -- LAPACK is a software package provided by Univ. of Tennessee, --
337 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
338 * June 2016
339 *
340 * .. Scalar Arguments ..
341  CHARACTER jobu1, jobu2, jobv1t, jobv2t, trans
342  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, lrwork, m, p, q
343 * ..
344 * .. Array Arguments ..
345  DOUBLE PRECISION b11d( * ), b11e( * ), b12d( * ), b12e( * ),
346  $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
347  $ phi( * ), theta( * ), rwork( * )
348  COMPLEX*16 u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
349  $ v2t( ldv2t, * )
350 * ..
351 *
352 * ===================================================================
353 *
354 * .. Parameters ..
355  INTEGER maxitr
356  parameter ( maxitr = 6 )
357  DOUBLE PRECISION hundred, meighth, one, piover2, ten, zero
358  parameter ( hundred = 100.0d0, meighth = -0.125d0,
359  $ one = 1.0d0, piover2 = 1.57079632679489662d0,
360  $ ten = 10.0d0, zero = 0.0d0 )
361  COMPLEX*16 negonecomplex
362  parameter ( negonecomplex = (-1.0d0,0.0d0) )
363 * ..
364 * .. Local Scalars ..
365  LOGICAL colmajor, lquery, restart11, restart12,
366  $ restart21, restart22, wantu1, wantu2, wantv1t,
367  $ wantv2t
368  INTEGER i, imin, imax, iter, iu1cs, iu1sn, iu2cs,
369  $ iu2sn, iv1tcs, iv1tsn, iv2tcs, iv2tsn, j,
370  $ lrworkmin, lrworkopt, maxit, mini
371  DOUBLE PRECISION b11bulge, b12bulge, b21bulge, b22bulge, dummy,
372  $ eps, mu, nu, r, sigma11, sigma21,
373  $ temp, thetamax, thetamin, thresh, tol, tolmul,
374  $ unfl, x1, x2, y1, y2
375 *
376  EXTERNAL dlartgp, dlartgs, dlas2, xerbla, zlasr, zscal,
377  $ zswap
378 * ..
379 * .. External Functions ..
380  DOUBLE PRECISION dlamch
381  LOGICAL lsame
382  EXTERNAL lsame, dlamch
383 * ..
384 * .. Intrinsic Functions ..
385  INTRINSIC abs, atan2, cos, max, min, sin, sqrt
386 * ..
387 * .. Executable Statements ..
388 *
389 * Test input arguments
390 *
391  info = 0
392  lquery = lrwork .EQ. -1
393  wantu1 = lsame( jobu1, 'Y' )
394  wantu2 = lsame( jobu2, 'Y' )
395  wantv1t = lsame( jobv1t, 'Y' )
396  wantv2t = lsame( jobv2t, 'Y' )
397  colmajor = .NOT. lsame( trans, 'T' )
398 *
399  IF( m .LT. 0 ) THEN
400  info = -6
401  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
402  info = -7
403  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
404  info = -8
405  ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q ) THEN
406  info = -8
407  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
408  info = -12
409  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
410  info = -14
411  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
412  info = -16
413  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
414  info = -18
415  END IF
416 *
417 * Quick return if Q = 0
418 *
419  IF( info .EQ. 0 .AND. q .EQ. 0 ) THEN
420  lrworkmin = 1
421  rwork(1) = lrworkmin
422  RETURN
423  END IF
424 *
425 * Compute workspace
426 *
427  IF( info .EQ. 0 ) THEN
428  iu1cs = 1
429  iu1sn = iu1cs + q
430  iu2cs = iu1sn + q
431  iu2sn = iu2cs + q
432  iv1tcs = iu2sn + q
433  iv1tsn = iv1tcs + q
434  iv2tcs = iv1tsn + q
435  iv2tsn = iv2tcs + q
436  lrworkopt = iv2tsn + q - 1
437  lrworkmin = lrworkopt
438  rwork(1) = lrworkopt
439  IF( lrwork .LT. lrworkmin .AND. .NOT. lquery ) THEN
440  info = -28
441  END IF
442  END IF
443 *
444  IF( info .NE. 0 ) THEN
445  CALL xerbla( 'ZBBCSD', -info )
446  RETURN
447  ELSE IF( lquery ) THEN
448  RETURN
449  END IF
450 *
451 * Get machine constants
452 *
453  eps = dlamch( 'Epsilon' )
454  unfl = dlamch( 'Safe minimum' )
455  tolmul = max( ten, min( hundred, eps**meighth ) )
456  tol = tolmul*eps
457  thresh = max( tol, maxitr*q*q*unfl )
458 *
459 * Test for negligible sines or cosines
460 *
461  DO i = 1, q
462  IF( theta(i) .LT. thresh ) THEN
463  theta(i) = zero
464  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
465  theta(i) = piover2
466  END IF
467  END DO
468  DO i = 1, q-1
469  IF( phi(i) .LT. thresh ) THEN
470  phi(i) = zero
471  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
472  phi(i) = piover2
473  END IF
474  END DO
475 *
476 * Initial deflation
477 *
478  imax = q
479  DO WHILE( imax .GT. 1 )
480  IF( phi(imax-1) .NE. zero ) THEN
481  EXIT
482  END IF
483  imax = imax - 1
484  END DO
485  imin = imax - 1
486  IF ( imin .GT. 1 ) THEN
487  DO WHILE( phi(imin-1) .NE. zero )
488  imin = imin - 1
489  IF ( imin .LE. 1 ) EXIT
490  END DO
491  END IF
492 *
493 * Initialize iteration counter
494 *
495  maxit = maxitr*q*q
496  iter = 0
497 *
498 * Begin main iteration loop
499 *
500  DO WHILE( imax .GT. 1 )
501 *
502 * Compute the matrix entries
503 *
504  b11d(imin) = cos( theta(imin) )
505  b21d(imin) = -sin( theta(imin) )
506  DO i = imin, imax - 1
507  b11e(i) = -sin( theta(i) ) * sin( phi(i) )
508  b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
509  b12d(i) = sin( theta(i) ) * cos( phi(i) )
510  b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
511  b21e(i) = -cos( theta(i) ) * sin( phi(i) )
512  b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
513  b22d(i) = cos( theta(i) ) * cos( phi(i) )
514  b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
515  END DO
516  b12d(imax) = sin( theta(imax) )
517  b22d(imax) = cos( theta(imax) )
518 *
519 * Abort if not converging; otherwise, increment ITER
520 *
521  IF( iter .GT. maxit ) THEN
522  info = 0
523  DO i = 1, q
524  IF( phi(i) .NE. zero )
525  $ info = info + 1
526  END DO
527  RETURN
528  END IF
529 *
530  iter = iter + imax - imin
531 *
532 * Compute shifts
533 *
534  thetamax = theta(imin)
535  thetamin = theta(imin)
536  DO i = imin+1, imax
537  IF( theta(i) > thetamax )
538  $ thetamax = theta(i)
539  IF( theta(i) < thetamin )
540  $ thetamin = theta(i)
541  END DO
542 *
543  IF( thetamax .GT. piover2 - thresh ) THEN
544 *
545 * Zero on diagonals of B11 and B22; induce deflation with a
546 * zero shift
547 *
548  mu = zero
549  nu = one
550 *
551  ELSE IF( thetamin .LT. thresh ) THEN
552 *
553 * Zero on diagonals of B12 and B22; induce deflation with a
554 * zero shift
555 *
556  mu = one
557  nu = zero
558 *
559  ELSE
560 *
561 * Compute shifts for B11 and B21 and use the lesser
562 *
563  CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
564  $ dummy )
565  CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), 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), rwork(iv2tsn+i-1-1),
723  $ rwork(iv2tcs+i-1-1), r )
724  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
725  CALL dlartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
726  $ rwork(iv2tcs+i-1-1), r )
727  ELSE IF( nu .LT. mu ) THEN
728  CALL dlartgs( b12e(i-1), b12d(i), nu,
729  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
730  ELSE
731  CALL dlartgs( b22e(i-1), b22d(i), mu,
732  $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
733  END IF
734 *
735  temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
736  b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
737  $ rwork(iv1tsn+i-1)*b11d(i)
738  b11d(i) = temp
739  b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
740  b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
741  temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
742  b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
743  $ rwork(iv1tsn+i-1)*b21d(i)
744  b21d(i) = temp
745  b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
746  b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
747  temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
748  $ rwork(iv2tsn+i-1-1)*b12d(i)
749  b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
750  $ rwork(iv2tsn+i-1-1)*b12e(i-1)
751  b12e(i-1) = temp
752  b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
753  b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
754  temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
755  $ rwork(iv2tsn+i-1-1)*b22d(i)
756  b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
757  $ rwork(iv2tsn+i-1-1)*b22e(i-1)
758  b22e(i-1) = temp
759  b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
760  b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
761 *
762 * Compute THETA(I)
763 *
764  x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
765  x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
766  y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
767  y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
768 *
769  theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
770 *
771 * Determine if there are bulges to chase or if a new direct
772 * summand has been reached
773 *
774  restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
775  restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
776  restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
777  restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
778 *
779 * If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
780 * B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
781 * chasing by applying the original shift again.
782 *
783  IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
784  CALL dlartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
785  $ r )
786  ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
787  CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
788  $ rwork(iu1cs+i-1), r )
789  ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
790  CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
791  $ rwork(iu1cs+i-1), r )
792  ELSE IF( mu .LE. nu ) THEN
793  CALL dlartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
794  $ rwork(iu1sn+i-1) )
795  ELSE
796  CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
797  $ rwork(iu1sn+i-1) )
798  END IF
799  IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
800  CALL dlartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
801  $ r )
802  ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
803  CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
804  $ rwork(iu2cs+i-1), r )
805  ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
806  CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
807  $ rwork(iu2cs+i-1), r )
808  ELSE IF( nu .LT. mu ) THEN
809  CALL dlartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
810  $ rwork(iu2sn+i-1) )
811  ELSE
812  CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
813  $ rwork(iu2sn+i-1) )
814  END IF
815  rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
816  rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
817 *
818  temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
819  b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
820  $ rwork(iu1sn+i-1)*b11e(i)
821  b11e(i) = temp
822  IF( i .LT. imax - 1 ) THEN
823  b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
824  b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
825  END IF
826  temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
827  b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
828  $ rwork(iu2sn+i-1)*b21e(i)
829  b21e(i) = temp
830  IF( i .LT. imax - 1 ) THEN
831  b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
832  b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
833  END IF
834  temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
835  b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
836  $ rwork(iu1sn+i-1)*b12d(i)
837  b12d(i) = temp
838  b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
839  b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
840  temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
841  b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
842  $ rwork(iu2sn+i-1)*b22d(i)
843  b22d(i) = temp
844  b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
845  b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
846 *
847  END DO
848 *
849 * Compute PHI(IMAX-1)
850 *
851  x1 = sin(theta(imax-1))*b11e(imax-1) +
852  $ cos(theta(imax-1))*b21e(imax-1)
853  y1 = sin(theta(imax-1))*b12d(imax-1) +
854  $ cos(theta(imax-1))*b22d(imax-1)
855  y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
856 *
857  phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
858 *
859 * Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
860 *
861  restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
862  restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
863 *
864  IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
865  CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
866  $ rwork(iv2tcs+imax-1-1), r )
867  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
868  CALL dlartgp( b12bulge, b12d(imax-1),
869  $ rwork(iv2tsn+imax-1-1),
870  $ rwork(iv2tcs+imax-1-1), r )
871  ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
872  CALL dlartgp( b22bulge, b22d(imax-1),
873  $ rwork(iv2tsn+imax-1-1),
874  $ rwork(iv2tcs+imax-1-1), r )
875  ELSE IF( nu .LT. mu ) THEN
876  CALL dlartgs( b12e(imax-1), b12d(imax), nu,
877  $ rwork(iv2tcs+imax-1-1),
878  $ rwork(iv2tsn+imax-1-1) )
879  ELSE
880  CALL dlartgs( b22e(imax-1), b22d(imax), mu,
881  $ rwork(iv2tcs+imax-1-1),
882  $ rwork(iv2tsn+imax-1-1) )
883  END IF
884 *
885  temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
886  $ rwork(iv2tsn+imax-1-1)*b12d(imax)
887  b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
888  $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
889  b12e(imax-1) = temp
890  temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
891  $ rwork(iv2tsn+imax-1-1)*b22d(imax)
892  b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
893  $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
894  b22e(imax-1) = temp
895 *
896 * Update singular vectors
897 *
898  IF( wantu1 ) THEN
899  IF( colmajor ) THEN
900  CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
901  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
902  $ u1(1,imin), ldu1 )
903  ELSE
904  CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
905  $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
906  $ u1(imin,1), ldu1 )
907  END IF
908  END IF
909  IF( wantu2 ) THEN
910  IF( colmajor ) THEN
911  CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
912  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
913  $ u2(1,imin), ldu2 )
914  ELSE
915  CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
916  $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
917  $ u2(imin,1), ldu2 )
918  END IF
919  END IF
920  IF( wantv1t ) THEN
921  IF( colmajor ) THEN
922  CALL zlasr( 'L', 'V', 'F', imax-imin+1, q,
923  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
924  $ v1t(imin,1), ldv1t )
925  ELSE
926  CALL zlasr( 'R', 'V', 'F', q, imax-imin+1,
927  $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
928  $ v1t(1,imin), ldv1t )
929  END IF
930  END IF
931  IF( wantv2t ) THEN
932  IF( colmajor ) THEN
933  CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-q,
934  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
935  $ v2t(imin,1), ldv2t )
936  ELSE
937  CALL zlasr( 'R', 'V', 'F', m-q, imax-imin+1,
938  $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
939  $ v2t(1,imin), ldv2t )
940  END IF
941  END IF
942 *
943 * Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
944 *
945  IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
946  b11d(imax) = -b11d(imax)
947  b21d(imax) = -b21d(imax)
948  IF( wantv1t ) THEN
949  IF( colmajor ) THEN
950  CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
951  ELSE
952  CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
953  END IF
954  END IF
955  END IF
956 *
957 * Compute THETA(IMAX)
958 *
959  x1 = cos(phi(imax-1))*b11d(imax) +
960  $ sin(phi(imax-1))*b12e(imax-1)
961  y1 = cos(phi(imax-1))*b21d(imax) +
962  $ sin(phi(imax-1))*b22e(imax-1)
963 *
964  theta(imax) = atan2( abs(y1), abs(x1) )
965 *
966 * Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
967 * and B22(IMAX,IMAX-1)
968 *
969  IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
970  b12d(imax) = -b12d(imax)
971  IF( wantu1 ) THEN
972  IF( colmajor ) THEN
973  CALL zscal( p, negonecomplex, u1(1,imax), 1 )
974  ELSE
975  CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
976  END IF
977  END IF
978  END IF
979  IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
980  b22d(imax) = -b22d(imax)
981  IF( wantu2 ) THEN
982  IF( colmajor ) THEN
983  CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
984  ELSE
985  CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
986  END IF
987  END IF
988  END IF
989 *
990 * Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
991 *
992  IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
993  IF( wantv2t ) THEN
994  IF( colmajor ) THEN
995  CALL zscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
996  ELSE
997  CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
998  END IF
999  END IF
1000  END IF
1001 *
1002 * Test for negligible sines or cosines
1003 *
1004  DO i = imin, imax
1005  IF( theta(i) .LT. thresh ) THEN
1006  theta(i) = zero
1007  ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1008  theta(i) = piover2
1009  END IF
1010  END DO
1011  DO i = imin, imax-1
1012  IF( phi(i) .LT. thresh ) THEN
1013  phi(i) = zero
1014  ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1015  phi(i) = piover2
1016  END IF
1017  END DO
1018 *
1019 * Deflate
1020 *
1021  IF (imax .GT. 1) THEN
1022  DO WHILE( phi(imax-1) .EQ. zero )
1023  imax = imax - 1
1024  IF (imax .LE. 1) EXIT
1025  END DO
1026  END IF
1027  IF( imin .GT. imax - 1 )
1028  $ imin = imax - 1
1029  IF (imin .GT. 1) THEN
1030  DO WHILE (phi(imin-1) .NE. zero)
1031  imin = imin - 1
1032  IF (imin .LE. 1) EXIT
1033  END DO
1034  END IF
1035 *
1036 * Repeat main iteration loop
1037 *
1038  END DO
1039 *
1040 * Postprocessing: order THETA from least to greatest
1041 *
1042  DO i = 1, q
1043 *
1044  mini = i
1045  thetamin = theta(i)
1046  DO j = i+1, q
1047  IF( theta(j) .LT. thetamin ) THEN
1048  mini = j
1049  thetamin = theta(j)
1050  END IF
1051  END DO
1052 *
1053  IF( mini .NE. i ) THEN
1054  theta(mini) = theta(i)
1055  theta(i) = thetamin
1056  IF( colmajor ) THEN
1057  IF( wantu1 )
1058  $ CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1059  IF( wantu2 )
1060  $ CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1061  IF( wantv1t )
1062  $ CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1063  IF( wantv2t )
1064  $ CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1065  $ ldv2t )
1066  ELSE
1067  IF( wantu1 )
1068  $ CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1069  IF( wantu2 )
1070  $ CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1071  IF( wantv1t )
1072  $ CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1073  IF( wantv2t )
1074  $ CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1075  END IF
1076  END IF
1077 *
1078  END DO
1079 *
1080  RETURN
1081 *
1082 * End of ZBBCSD
1083 *
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:202
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:52
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:92
subroutine dlartgp(F, G, CS, SN, R)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: dlartgp.f:97
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: dlas2.f:109
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:54

Here is the call graph for this function:

Here is the caller graph for this function: