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

◆ dbbcsd()

subroutine dbbcsd ( character  jobu1,
character  jobu2,
character  jobv1t,
character  jobv2t,
character  trans,
integer  m,
integer  p,
integer  q,
double precision, dimension( * )  theta,
double precision, dimension( * )  phi,
double precision, dimension( ldu1, * )  u1,
integer  ldu1,
double precision, dimension( ldu2, * )  u2,
integer  ldu2,
double precision, dimension( ldv1t, * )  v1t,
integer  ldv1t,
double precision, 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( * )  work,
integer  lwork,
integer  info 
)

DBBCSD

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

Purpose:
 DBBCSD computes the CS decomposition of an orthogonal 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 |    ]**T
                 = [---------] [---------------] [---------]   .
                   [    | 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 DORCSD for details.)

 The bidiagonal matrices B11, B12, B21, and B22 are represented
 implicitly by angles THETA(1:Q) and PHI(1:Q-1).

 The orthogonal 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 orthogonal 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDV1T,Q)
          On entry, a Q-by-Q matrix. On exit, V1T is premultiplied
          by the 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 DOUBLE PRECISION array, dimension (LDV2T,M-Q)
          On entry, an (M-Q)-by-(M-Q) matrix. On exit, V2T is
          premultiplied by the 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 DBBCSD converges, B11D contains the cosines of THETA(1),
          ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B11E contains zeros. If DBBCSD 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 DBBCSD converges, B12D contains the negative sines of
          THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B12E contains zeros. If DBBCSD 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 DBBCSD converges, B21D contains the negative sines of
          THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B21E contains zeros. If DBBCSD 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 DBBCSD converges, B22D contains the negative sines of
          THETA(1), ..., THETA(Q). If DBBCSD 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 DBBCSD converges, B22E contains zeros. If DBBCSD fails
          to converge, then B22E contains the subdiagonal of the
          partially reduced bottom-right block.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= MAX(1,8*Q).

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