LAPACK 3.12.0
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 328 of file zbbcsd.f.

332*
333* -- LAPACK computational routine --
334* -- LAPACK is a software package provided by Univ. of Tennessee, --
335* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336*
337* .. Scalar Arguments ..
338 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
339 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LRWORK, M, P, Q
340* ..
341* .. Array Arguments ..
342 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344 $ PHI( * ), THETA( * ), RWORK( * )
345 COMPLEX*16 U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
346 $ V2T( LDV2T, * )
347* ..
348*
349* ===================================================================
350*
351* .. Parameters ..
352 INTEGER MAXITR
353 parameter( maxitr = 6 )
354 DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, TEN, ZERO
355 parameter( hundred = 100.0d0, meighth = -0.125d0,
356 $ one = 1.0d0, ten = 10.0d0, zero = 0.0d0 )
357 COMPLEX*16 NEGONECOMPLEX
358 parameter( negonecomplex = (-1.0d0,0.0d0) )
359 DOUBLE PRECISION PIOVER2
360 parameter( piover2 = 1.57079632679489661923132169163975144210d0 )
361* ..
362* .. Local Scalars ..
363 LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
364 $ RESTART21, RESTART22, WANTU1, WANTU2, WANTV1T,
365 $ WANTV2T
366 INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
367 $ IU2SN, IV1TCS, IV1TSN, IV2TCS, IV2TSN, J,
368 $ LRWORKMIN, LRWORKOPT, MAXIT, MINI
369 DOUBLE PRECISION B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
370 $ EPS, MU, NU, R, SIGMA11, SIGMA21,
371 $ TEMP, THETAMAX, THETAMIN, THRESH, TOL, TOLMUL,
372 $ UNFL, X1, X2, Y1, Y2
373*
374 EXTERNAL dlartgp, dlartgs, dlas2, xerbla, zlasr, 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), sigma11,
562 $ dummy )
563 CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
564 $ dummy )
565*
566 IF( sigma11 .LE. sigma21 ) THEN
567 mu = sigma11
568 nu = sqrt( one - mu**2 )
569 IF( mu .LT. thresh ) THEN
570 mu = zero
571 nu = one
572 END IF
573 ELSE
574 nu = sigma21
575 mu = sqrt( 1.0 - nu**2 )
576 IF( nu .LT. thresh ) THEN
577 mu = one
578 nu = zero
579 END IF
580 END IF
581 END IF
582*
583* Rotate to produce bulges in B11 and B21
584*
585 IF( mu .LE. nu ) THEN
586 CALL dlartgs( b11d(imin), b11e(imin), mu,
587 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
588 ELSE
589 CALL dlartgs( b21d(imin), b21e(imin), nu,
590 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
591 END IF
592*
593 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
594 $ rwork(iv1tsn+imin-1)*b11e(imin)
595 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
596 $ rwork(iv1tsn+imin-1)*b11d(imin)
597 b11d(imin) = temp
598 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
599 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
600 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
601 $ rwork(iv1tsn+imin-1)*b21e(imin)
602 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
603 $ rwork(iv1tsn+imin-1)*b21d(imin)
604 b21d(imin) = temp
605 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
606 b21d(imin+1) = rwork(iv1tcs+imin-1)*b21d(imin+1)
607*
608* Compute THETA(IMIN)
609*
610 theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
611 $ sqrt( b11d(imin)**2+b11bulge**2 ) )
612*
613* Chase the bulges in B11(IMIN+1,IMIN) and B21(IMIN+1,IMIN)
614*
615 IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 ) THEN
616 CALL dlartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
617 $ rwork(iu1cs+imin-1), r )
618 ELSE IF( mu .LE. nu ) THEN
619 CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
620 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
621 ELSE
622 CALL dlartgs( b12d( imin ), b12e( imin ), nu,
623 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
624 END IF
625 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
626 CALL dlartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
627 $ rwork(iu2cs+imin-1), r )
628 ELSE IF( nu .LT. mu ) THEN
629 CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
630 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
631 ELSE
632 CALL dlartgs( b22d(imin), b22e(imin), mu,
633 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
634 END IF
635 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
636 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
637*
638 temp = rwork(iu1cs+imin-1)*b11e(imin) +
639 $ rwork(iu1sn+imin-1)*b11d(imin+1)
640 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
641 $ rwork(iu1sn+imin-1)*b11e(imin)
642 b11e(imin) = temp
643 IF( imax .GT. imin+1 ) THEN
644 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
645 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
646 END IF
647 temp = rwork(iu1cs+imin-1)*b12d(imin) +
648 $ rwork(iu1sn+imin-1)*b12e(imin)
649 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
650 $ rwork(iu1sn+imin-1)*b12d(imin)
651 b12d(imin) = temp
652 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
653 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
654 temp = rwork(iu2cs+imin-1)*b21e(imin) +
655 $ rwork(iu2sn+imin-1)*b21d(imin+1)
656 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
657 $ rwork(iu2sn+imin-1)*b21e(imin)
658 b21e(imin) = temp
659 IF( imax .GT. imin+1 ) THEN
660 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
661 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
662 END IF
663 temp = rwork(iu2cs+imin-1)*b22d(imin) +
664 $ rwork(iu2sn+imin-1)*b22e(imin)
665 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
666 $ rwork(iu2sn+imin-1)*b22d(imin)
667 b22d(imin) = temp
668 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
669 b22d(imin+1) = rwork(iu2cs+imin-1)*b22d(imin+1)
670*
671* Inner loop: chase bulges from B11(IMIN,IMIN+2),
672* B12(IMIN,IMIN+1), B21(IMIN,IMIN+2), and B22(IMIN,IMIN+1) to
673* bottom-right
674*
675 DO i = imin+1, imax-1
676*
677* Compute PHI(I-1)
678*
679 x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
680 x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
681 y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
682 y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
683*
684 phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
685*
686* Determine if there are bulges to chase or if a new direct
687* summand has been reached
688*
689 restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
690 restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
691 restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
692 restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
693*
694* If possible, chase bulges from B11(I-1,I+1), B12(I-1,I),
695* B21(I-1,I+1), and B22(I-1,I). If necessary, restart bulge-
696* chasing by applying the original shift again.
697*
698 IF( .NOT. restart11 .AND. .NOT. restart21 ) THEN
699 CALL dlartgp( x2, x1, rwork(iv1tsn+i-1),
700 $ rwork(iv1tcs+i-1), r )
701 ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
702 CALL dlartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
703 $ rwork(iv1tcs+i-1), r )
704 ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
705 CALL dlartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
706 $ rwork(iv1tcs+i-1), r )
707 ELSE IF( mu .LE. nu ) THEN
708 CALL dlartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
709 $ rwork(iv1tsn+i-1) )
710 ELSE
711 CALL dlartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
712 $ rwork(iv1tsn+i-1) )
713 END IF
714 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
715 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
716 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
717 CALL dlartgp( y2, y1, rwork(iv2tsn+i-1-1),
718 $ rwork(iv2tcs+i-1-1), r )
719 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
720 CALL dlartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
721 $ rwork(iv2tcs+i-1-1), r )
722 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
723 CALL dlartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
724 $ rwork(iv2tcs+i-1-1), r )
725 ELSE IF( nu .LT. mu ) THEN
726 CALL dlartgs( b12e(i-1), b12d(i), nu,
727 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
728 ELSE
729 CALL dlartgs( b22e(i-1), b22d(i), mu,
730 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
731 END IF
732*
733 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
734 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
735 $ rwork(iv1tsn+i-1)*b11d(i)
736 b11d(i) = temp
737 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
738 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
739 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
740 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
741 $ rwork(iv1tsn+i-1)*b21d(i)
742 b21d(i) = temp
743 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
744 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
745 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
746 $ rwork(iv2tsn+i-1-1)*b12d(i)
747 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
748 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
749 b12e(i-1) = temp
750 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
751 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
752 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
753 $ rwork(iv2tsn+i-1-1)*b22d(i)
754 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
755 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
756 b22e(i-1) = temp
757 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
758 b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
759*
760* Compute THETA(I)
761*
762 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
763 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
764 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
765 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
766*
767 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
768*
769* Determine if there are bulges to chase or if a new direct
770* summand has been reached
771*
772 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
773 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
774 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
775 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
776*
777* If possible, chase bulges from B11(I+1,I), B12(I+1,I-1),
778* B21(I+1,I), and B22(I+1,I-1). If necessary, restart bulge-
779* chasing by applying the original shift again.
780*
781 IF( .NOT. restart11 .AND. .NOT. restart12 ) THEN
782 CALL dlartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
783 $ r )
784 ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
785 CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
786 $ rwork(iu1cs+i-1), r )
787 ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
788 CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
789 $ rwork(iu1cs+i-1), r )
790 ELSE IF( mu .LE. nu ) THEN
791 CALL dlartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
792 $ rwork(iu1sn+i-1) )
793 ELSE
794 CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
795 $ rwork(iu1sn+i-1) )
796 END IF
797 IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
798 CALL dlartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
799 $ r )
800 ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
801 CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
802 $ rwork(iu2cs+i-1), r )
803 ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
804 CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
805 $ rwork(iu2cs+i-1), r )
806 ELSE IF( nu .LT. mu ) THEN
807 CALL dlartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
808 $ rwork(iu2sn+i-1) )
809 ELSE
810 CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
811 $ rwork(iu2sn+i-1) )
812 END IF
813 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
814 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
815*
816 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
817 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
818 $ rwork(iu1sn+i-1)*b11e(i)
819 b11e(i) = temp
820 IF( i .LT. imax - 1 ) THEN
821 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
822 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
823 END IF
824 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
825 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
826 $ rwork(iu2sn+i-1)*b21e(i)
827 b21e(i) = temp
828 IF( i .LT. imax - 1 ) THEN
829 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
830 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
831 END IF
832 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
833 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
834 $ rwork(iu1sn+i-1)*b12d(i)
835 b12d(i) = temp
836 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
837 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
838 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
839 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
840 $ rwork(iu2sn+i-1)*b22d(i)
841 b22d(i) = temp
842 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
843 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
844*
845 END DO
846*
847* Compute PHI(IMAX-1)
848*
849 x1 = sin(theta(imax-1))*b11e(imax-1) +
850 $ cos(theta(imax-1))*b21e(imax-1)
851 y1 = sin(theta(imax-1))*b12d(imax-1) +
852 $ cos(theta(imax-1))*b22d(imax-1)
853 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
854*
855 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
856*
857* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
858*
859 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
860 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
861*
862 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
863 CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
864 $ rwork(iv2tcs+imax-1-1), r )
865 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
866 CALL dlartgp( b12bulge, b12d(imax-1),
867 $ rwork(iv2tsn+imax-1-1),
868 $ rwork(iv2tcs+imax-1-1), r )
869 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
870 CALL dlartgp( b22bulge, b22d(imax-1),
871 $ rwork(iv2tsn+imax-1-1),
872 $ rwork(iv2tcs+imax-1-1), r )
873 ELSE IF( nu .LT. mu ) THEN
874 CALL dlartgs( b12e(imax-1), b12d(imax), nu,
875 $ rwork(iv2tcs+imax-1-1),
876 $ rwork(iv2tsn+imax-1-1) )
877 ELSE
878 CALL dlartgs( b22e(imax-1), b22d(imax), mu,
879 $ rwork(iv2tcs+imax-1-1),
880 $ rwork(iv2tsn+imax-1-1) )
881 END IF
882*
883 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
884 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
885 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
886 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
887 b12e(imax-1) = temp
888 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
889 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
890 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
891 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
892 b22e(imax-1) = temp
893*
894* Update singular vectors
895*
896 IF( wantu1 ) THEN
897 IF( colmajor ) THEN
898 CALL zlasr( 'R', 'V', 'F', p, imax-imin+1,
899 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
900 $ u1(1,imin), ldu1 )
901 ELSE
902 CALL zlasr( 'L', 'V', 'F', imax-imin+1, p,
903 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
904 $ u1(imin,1), ldu1 )
905 END IF
906 END IF
907 IF( wantu2 ) THEN
908 IF( colmajor ) THEN
909 CALL zlasr( 'R', 'V', 'F', m-p, imax-imin+1,
910 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
911 $ u2(1,imin), ldu2 )
912 ELSE
913 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-p,
914 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
915 $ u2(imin,1), ldu2 )
916 END IF
917 END IF
918 IF( wantv1t ) THEN
919 IF( colmajor ) THEN
920 CALL zlasr( 'L', 'V', 'F', imax-imin+1, q,
921 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
922 $ v1t(imin,1), ldv1t )
923 ELSE
924 CALL zlasr( 'R', 'V', 'F', q, imax-imin+1,
925 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
926 $ v1t(1,imin), ldv1t )
927 END IF
928 END IF
929 IF( wantv2t ) THEN
930 IF( colmajor ) THEN
931 CALL zlasr( 'L', 'V', 'F', imax-imin+1, m-q,
932 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
933 $ v2t(imin,1), ldv2t )
934 ELSE
935 CALL zlasr( 'R', 'V', 'F', m-q, imax-imin+1,
936 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
937 $ v2t(1,imin), ldv2t )
938 END IF
939 END IF
940*
941* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
942*
943 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
944 b11d(imax) = -b11d(imax)
945 b21d(imax) = -b21d(imax)
946 IF( wantv1t ) THEN
947 IF( colmajor ) THEN
948 CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
949 ELSE
950 CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
951 END IF
952 END IF
953 END IF
954*
955* Compute THETA(IMAX)
956*
957 x1 = cos(phi(imax-1))*b11d(imax) +
958 $ sin(phi(imax-1))*b12e(imax-1)
959 y1 = cos(phi(imax-1))*b21d(imax) +
960 $ sin(phi(imax-1))*b22e(imax-1)
961*
962 theta(imax) = atan2( abs(y1), abs(x1) )
963*
964* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
965* and B22(IMAX,IMAX-1)
966*
967 IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
968 b12d(imax) = -b12d(imax)
969 IF( wantu1 ) THEN
970 IF( colmajor ) THEN
971 CALL zscal( p, negonecomplex, u1(1,imax), 1 )
972 ELSE
973 CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
974 END IF
975 END IF
976 END IF
977 IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
978 b22d(imax) = -b22d(imax)
979 IF( wantu2 ) THEN
980 IF( colmajor ) THEN
981 CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
982 ELSE
983 CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
984 END IF
985 END IF
986 END IF
987*
988* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
989*
990 IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
991 IF( wantv2t ) THEN
992 IF( colmajor ) THEN
993 CALL zscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
994 ELSE
995 CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
996 END IF
997 END IF
998 END IF
999*
1000* Test for negligible sines or cosines
1001*
1002 DO i = imin, imax
1003 IF( theta(i) .LT. thresh ) THEN
1004 theta(i) = zero
1005 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1006 theta(i) = piover2
1007 END IF
1008 END DO
1009 DO i = imin, imax-1
1010 IF( phi(i) .LT. thresh ) THEN
1011 phi(i) = zero
1012 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1013 phi(i) = piover2
1014 END IF
1015 END DO
1016*
1017* Deflate
1018*
1019 IF (imax .GT. 1) THEN
1020 DO WHILE( phi(imax-1) .EQ. zero )
1021 imax = imax - 1
1022 IF (imax .LE. 1) EXIT
1023 END DO
1024 END IF
1025 IF( imin .GT. imax - 1 )
1026 $ imin = imax - 1
1027 IF (imin .GT. 1) THEN
1028 DO WHILE (phi(imin-1) .NE. zero)
1029 imin = imin - 1
1030 IF (imin .LE. 1) EXIT
1031 END DO
1032 END IF
1033*
1034* Repeat main iteration loop
1035*
1036 END DO
1037*
1038* Postprocessing: order THETA from least to greatest
1039*
1040 DO i = 1, q
1041*
1042 mini = i
1043 thetamin = theta(i)
1044 DO j = i+1, q
1045 IF( theta(j) .LT. thetamin ) THEN
1046 mini = j
1047 thetamin = theta(j)
1048 END IF
1049 END DO
1050*
1051 IF( mini .NE. i ) THEN
1052 theta(mini) = theta(i)
1053 theta(i) = thetamin
1054 IF( colmajor ) THEN
1055 IF( wantu1 )
1056 $ CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1057 IF( wantu2 )
1058 $ CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1059 IF( wantv1t )
1060 $ CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1061 IF( wantv2t )
1062 $ CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1063 $ ldv2t )
1064 ELSE
1065 IF( wantu1 )
1066 $ CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1067 IF( wantu2 )
1068 $ CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1069 IF( wantv1t )
1070 $ CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1071 IF( wantv2t )
1072 $ CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1073 END IF
1074 END IF
1075*
1076 END DO
1077*
1078 RETURN
1079*
1080* End of ZBBCSD
1081*
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:95
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:90
dlas2
subroutine dlas2(f, g, h, ssmin, ssmax)
DLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition dlas2.f:105
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:200
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: