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

◆ cbbcsd()

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

CBBCSD

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

Purpose:
 CBBCSD 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 CUNCSD 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 REAL 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 REAL 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 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 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 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 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 REAL array, dimension (Q)
          When CBBCSD converges, B11D contains the cosines of THETA(1),
          ..., THETA(Q). If CBBCSD fails to converge, then B11D
          contains the diagonal of the partially reduced top-left
          block.
[out]B11E
          B11E is REAL array, dimension (Q-1)
          When CBBCSD converges, B11E contains zeros. If CBBCSD fails
          to converge, then B11E contains the superdiagonal of the
          partially reduced top-left block.
[out]B12D
          B12D is REAL array, dimension (Q)
          When CBBCSD converges, B12D contains the negative sines of
          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
          B12D contains the diagonal of the partially reduced top-right
          block.
[out]B12E
          B12E is REAL array, dimension (Q-1)
          When CBBCSD converges, B12E contains zeros. If CBBCSD fails
          to converge, then B12E contains the subdiagonal of the
          partially reduced top-right block.
[out]B21D
          B21D is REAL array, dimension (Q)
          When CBBCSD converges, B21D contains the negative sines of
          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
          B21D contains the diagonal of the partially reduced bottom-left
          block.
[out]B21E
          B21E is REAL array, dimension (Q-1)
          When CBBCSD converges, B21E contains zeros. If CBBCSD fails
          to converge, then B21E contains the subdiagonal of the
          partially reduced bottom-left block.
[out]B22D
          B22D is REAL array, dimension (Q)
          When CBBCSD converges, B22D contains the negative sines of
          THETA(1), ..., THETA(Q). If CBBCSD fails to converge, then
          B22D contains the diagonal of the partially reduced bottom-right
          block.
[out]B22E
          B22E is REAL array, dimension (Q-1)
          When CBBCSD converges, B22E contains zeros. If CBBCSD fails
          to converge, then B22E contains the subdiagonal of the
          partially reduced bottom-right block.
[out]RWORK
          RWORK is REAL 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 CBBCSD did not converge, INFO specifies the number
                of nonzero entries in PHI, and B11D, B11E, etc.,
                contain the partially reduced matrix.
Internal Parameters:
  TOLMUL  REAL, 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 cbbcsd.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 REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343 $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344 $ PHI( * ), THETA( * ), RWORK( * )
345 COMPLEX U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
346 $ V2T( LDV2T, * )
347* ..
348*
349* ===================================================================
350*
351* .. Parameters ..
352 INTEGER MAXITR
353 parameter( maxitr = 6 )
354 REAL HUNDRED, MEIGHTH, ONE, TEN, ZERO
355 parameter( hundred = 100.0e0, meighth = -0.125e0,
356 $ one = 1.0e0, ten = 10.0e0, zero = 0.0e0 )
357 COMPLEX NEGONECOMPLEX
358 parameter( negonecomplex = (-1.0e0,0.0e0) )
359 REAL PIOVER2
360 parameter( piover2 = 1.57079632679489661923132169163975144210e0 )
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 REAL 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 clasr, cscal, cswap, slartgp, slartgs, slas2,
376 $ xerbla
377* ..
378* .. External Functions ..
379 REAL SLAMCH
380 LOGICAL LSAME
381 EXTERNAL lsame, slamch
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 = lrwork .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 lrworkmin = 1
420 rwork(1) = lrworkmin
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 lrworkopt = iv2tsn + q - 1
436 lrworkmin = lrworkopt
437 rwork(1) = lrworkopt
438 IF( lrwork .LT. lrworkmin .AND. .NOT. lquery ) THEN
439 info = -28
440 END IF
441 END IF
442*
443 IF( info .NE. 0 ) THEN
444 CALL xerbla( 'CBBCSD', -info )
445 RETURN
446 ELSE IF( lquery ) THEN
447 RETURN
448 END IF
449*
450* Get machine constants
451*
452 eps = slamch( 'Epsilon' )
453 unfl = slamch( '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 slas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
563 $ dummy )
564 CALL slas2( 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 slartgs( b11d(imin), b11e(imin), mu,
588 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
589 ELSE
590 CALL slartgs( b21d(imin), b21e(imin), nu,
591 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
592 END IF
593*
594 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
595 $ rwork(iv1tsn+imin-1)*b11e(imin)
596 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
597 $ rwork(iv1tsn+imin-1)*b11d(imin)
598 b11d(imin) = temp
599 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
600 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
601 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
602 $ rwork(iv1tsn+imin-1)*b21e(imin)
603 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
604 $ rwork(iv1tsn+imin-1)*b21d(imin)
605 b21d(imin) = temp
606 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
607 b21d(imin+1) = rwork(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 slartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
618 $ rwork(iu1cs+imin-1), r )
619 ELSE IF( mu .LE. nu ) THEN
620 CALL slartgs( b11e( imin ), b11d( imin + 1 ), mu,
621 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
622 ELSE
623 CALL slartgs( b12d( imin ), b12e( imin ), nu,
624 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
625 END IF
626 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 ) THEN
627 CALL slartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
628 $ rwork(iu2cs+imin-1), r )
629 ELSE IF( nu .LT. mu ) THEN
630 CALL slartgs( b21e( imin ), b21d( imin + 1 ), nu,
631 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
632 ELSE
633 CALL slartgs( b22d(imin), b22e(imin), mu,
634 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
635 END IF
636 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
637 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
638*
639 temp = rwork(iu1cs+imin-1)*b11e(imin) +
640 $ rwork(iu1sn+imin-1)*b11d(imin+1)
641 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
642 $ rwork(iu1sn+imin-1)*b11e(imin)
643 b11e(imin) = temp
644 IF( imax .GT. imin+1 ) THEN
645 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
646 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
647 END IF
648 temp = rwork(iu1cs+imin-1)*b12d(imin) +
649 $ rwork(iu1sn+imin-1)*b12e(imin)
650 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
651 $ rwork(iu1sn+imin-1)*b12d(imin)
652 b12d(imin) = temp
653 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
654 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
655 temp = rwork(iu2cs+imin-1)*b21e(imin) +
656 $ rwork(iu2sn+imin-1)*b21d(imin+1)
657 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
658 $ rwork(iu2sn+imin-1)*b21e(imin)
659 b21e(imin) = temp
660 IF( imax .GT. imin+1 ) THEN
661 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
662 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
663 END IF
664 temp = rwork(iu2cs+imin-1)*b22d(imin) +
665 $ rwork(iu2sn+imin-1)*b22e(imin)
666 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
667 $ rwork(iu2sn+imin-1)*b22d(imin)
668 b22d(imin) = temp
669 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
670 b22d(imin+1) = rwork(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 slartgp( x2, x1, rwork(iv1tsn+i-1),
701 $ rwork(iv1tcs+i-1), r )
702 ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
703 CALL slartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
704 $ rwork(iv1tcs+i-1), r )
705 ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
706 CALL slartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
707 $ rwork(iv1tcs+i-1), r )
708 ELSE IF( mu .LE. nu ) THEN
709 CALL slartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
710 $ rwork(iv1tsn+i-1) )
711 ELSE
712 CALL slartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
713 $ rwork(iv1tsn+i-1) )
714 END IF
715 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
716 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
717 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
718 CALL slartgp( y2, y1, rwork(iv2tsn+i-1-1),
719 $ rwork(iv2tcs+i-1-1), r )
720 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
721 CALL slartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
722 $ rwork(iv2tcs+i-1-1), r )
723 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
724 CALL slartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
725 $ rwork(iv2tcs+i-1-1), r )
726 ELSE IF( nu .LT. mu ) THEN
727 CALL slartgs( b12e(i-1), b12d(i), nu,
728 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
729 ELSE
730 CALL slartgs( b22e(i-1), b22d(i), mu,
731 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
732 END IF
733*
734 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
735 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
736 $ rwork(iv1tsn+i-1)*b11d(i)
737 b11d(i) = temp
738 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
739 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
740 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
741 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
742 $ rwork(iv1tsn+i-1)*b21d(i)
743 b21d(i) = temp
744 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
745 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
746 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
747 $ rwork(iv2tsn+i-1-1)*b12d(i)
748 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
749 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
750 b12e(i-1) = temp
751 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
752 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
753 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
754 $ rwork(iv2tsn+i-1-1)*b22d(i)
755 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
756 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
757 b22e(i-1) = temp
758 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
759 b22e(i) = rwork(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 slartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
784 $ r )
785 ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
786 CALL slartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
787 $ rwork(iu1cs+i-1), r )
788 ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
789 CALL slartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
790 $ rwork(iu1cs+i-1), r )
791 ELSE IF( mu .LE. nu ) THEN
792 CALL slartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
793 $ rwork(iu1sn+i-1) )
794 ELSE
795 CALL slartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
796 $ rwork(iu1sn+i-1) )
797 END IF
798 IF( .NOT. restart21 .AND. .NOT. restart22 ) THEN
799 CALL slartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
800 $ r )
801 ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
802 CALL slartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
803 $ rwork(iu2cs+i-1), r )
804 ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
805 CALL slartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
806 $ rwork(iu2cs+i-1), r )
807 ELSE IF( nu .LT. mu ) THEN
808 CALL slartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
809 $ rwork(iu2sn+i-1) )
810 ELSE
811 CALL slartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
812 $ rwork(iu2sn+i-1) )
813 END IF
814 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
815 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
816*
817 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
818 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
819 $ rwork(iu1sn+i-1)*b11e(i)
820 b11e(i) = temp
821 IF( i .LT. imax - 1 ) THEN
822 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
823 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
824 END IF
825 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
826 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
827 $ rwork(iu2sn+i-1)*b21e(i)
828 b21e(i) = temp
829 IF( i .LT. imax - 1 ) THEN
830 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
831 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
832 END IF
833 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
834 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
835 $ rwork(iu1sn+i-1)*b12d(i)
836 b12d(i) = temp
837 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
838 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
839 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
840 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
841 $ rwork(iu2sn+i-1)*b22d(i)
842 b22d(i) = temp
843 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
844 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
845*
846 END DO
847*
848* Compute PHI(IMAX-1)
849*
850 x1 = sin(theta(imax-1))*b11e(imax-1) +
851 $ cos(theta(imax-1))*b21e(imax-1)
852 y1 = sin(theta(imax-1))*b12d(imax-1) +
853 $ cos(theta(imax-1))*b22d(imax-1)
854 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
855*
856 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
857*
858* Chase bulges from B12(IMAX-1,IMAX) and B22(IMAX-1,IMAX)
859*
860 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
861 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
862*
863 IF( .NOT. restart12 .AND. .NOT. restart22 ) THEN
864 CALL slartgp( y2, y1, rwork(iv2tsn+imax-1-1),
865 $ rwork(iv2tcs+imax-1-1), r )
866 ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
867 CALL slartgp( b12bulge, b12d(imax-1),
868 $ rwork(iv2tsn+imax-1-1),
869 $ rwork(iv2tcs+imax-1-1), r )
870 ELSE IF( restart12 .AND. .NOT. restart22 ) THEN
871 CALL slartgp( b22bulge, b22d(imax-1),
872 $ rwork(iv2tsn+imax-1-1),
873 $ rwork(iv2tcs+imax-1-1), r )
874 ELSE IF( nu .LT. mu ) THEN
875 CALL slartgs( b12e(imax-1), b12d(imax), nu,
876 $ rwork(iv2tcs+imax-1-1),
877 $ rwork(iv2tsn+imax-1-1) )
878 ELSE
879 CALL slartgs( b22e(imax-1), b22d(imax), mu,
880 $ rwork(iv2tcs+imax-1-1),
881 $ rwork(iv2tsn+imax-1-1) )
882 END IF
883*
884 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
885 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
886 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
887 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
888 b12e(imax-1) = temp
889 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
890 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
891 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
892 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
893 b22e(imax-1) = temp
894*
895* Update singular vectors
896*
897 IF( wantu1 ) THEN
898 IF( colmajor ) THEN
899 CALL clasr( 'R', 'V', 'F', p, imax-imin+1,
900 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
901 $ u1(1,imin), ldu1 )
902 ELSE
903 CALL clasr( 'L', 'V', 'F', imax-imin+1, p,
904 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
905 $ u1(imin,1), ldu1 )
906 END IF
907 END IF
908 IF( wantu2 ) THEN
909 IF( colmajor ) THEN
910 CALL clasr( 'R', 'V', 'F', m-p, imax-imin+1,
911 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
912 $ u2(1,imin), ldu2 )
913 ELSE
914 CALL clasr( 'L', 'V', 'F', imax-imin+1, m-p,
915 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
916 $ u2(imin,1), ldu2 )
917 END IF
918 END IF
919 IF( wantv1t ) THEN
920 IF( colmajor ) THEN
921 CALL clasr( 'L', 'V', 'F', imax-imin+1, q,
922 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
923 $ v1t(imin,1), ldv1t )
924 ELSE
925 CALL clasr( 'R', 'V', 'F', q, imax-imin+1,
926 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
927 $ v1t(1,imin), ldv1t )
928 END IF
929 END IF
930 IF( wantv2t ) THEN
931 IF( colmajor ) THEN
932 CALL clasr( 'L', 'V', 'F', imax-imin+1, m-q,
933 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
934 $ v2t(imin,1), ldv2t )
935 ELSE
936 CALL clasr( 'R', 'V', 'F', m-q, imax-imin+1,
937 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
938 $ v2t(1,imin), ldv2t )
939 END IF
940 END IF
941*
942* Fix signs on B11(IMAX-1,IMAX) and B21(IMAX-1,IMAX)
943*
944 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 ) THEN
945 b11d(imax) = -b11d(imax)
946 b21d(imax) = -b21d(imax)
947 IF( wantv1t ) THEN
948 IF( colmajor ) THEN
949 CALL cscal( q, negonecomplex, v1t(imax,1), ldv1t )
950 ELSE
951 CALL cscal( q, negonecomplex, v1t(1,imax), 1 )
952 END IF
953 END IF
954 END IF
955*
956* Compute THETA(IMAX)
957*
958 x1 = cos(phi(imax-1))*b11d(imax) +
959 $ sin(phi(imax-1))*b12e(imax-1)
960 y1 = cos(phi(imax-1))*b21d(imax) +
961 $ sin(phi(imax-1))*b22e(imax-1)
962*
963 theta(imax) = atan2( abs(y1), abs(x1) )
964*
965* Fix signs on B11(IMAX,IMAX), B12(IMAX,IMAX-1), B21(IMAX,IMAX),
966* and B22(IMAX,IMAX-1)
967*
968 IF( b11d(imax)+b12e(imax-1) .LT. 0 ) THEN
969 b12d(imax) = -b12d(imax)
970 IF( wantu1 ) THEN
971 IF( colmajor ) THEN
972 CALL cscal( p, negonecomplex, u1(1,imax), 1 )
973 ELSE
974 CALL cscal( p, negonecomplex, u1(imax,1), ldu1 )
975 END IF
976 END IF
977 END IF
978 IF( b21d(imax)+b22e(imax-1) .GT. 0 ) THEN
979 b22d(imax) = -b22d(imax)
980 IF( wantu2 ) THEN
981 IF( colmajor ) THEN
982 CALL cscal( m-p, negonecomplex, u2(1,imax), 1 )
983 ELSE
984 CALL cscal( m-p, negonecomplex, u2(imax,1), ldu2 )
985 END IF
986 END IF
987 END IF
988*
989* Fix signs on B12(IMAX,IMAX) and B22(IMAX,IMAX)
990*
991 IF( b12d(imax)+b22d(imax) .LT. 0 ) THEN
992 IF( wantv2t ) THEN
993 IF( colmajor ) THEN
994 CALL cscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
995 ELSE
996 CALL cscal( m-q, negonecomplex, v2t(1,imax), 1 )
997 END IF
998 END IF
999 END IF
1000*
1001* Test for negligible sines or cosines
1002*
1003 DO i = imin, imax
1004 IF( theta(i) .LT. thresh ) THEN
1005 theta(i) = zero
1006 ELSE IF( theta(i) .GT. piover2-thresh ) THEN
1007 theta(i) = piover2
1008 END IF
1009 END DO
1010 DO i = imin, imax-1
1011 IF( phi(i) .LT. thresh ) THEN
1012 phi(i) = zero
1013 ELSE IF( phi(i) .GT. piover2-thresh ) THEN
1014 phi(i) = piover2
1015 END IF
1016 END DO
1017*
1018* Deflate
1019*
1020 IF (imax .GT. 1) THEN
1021 DO WHILE( phi(imax-1) .EQ. zero )
1022 imax = imax - 1
1023 IF (imax .LE. 1) EXIT
1024 END DO
1025 END IF
1026 IF( imin .GT. imax - 1 )
1027 $ imin = imax - 1
1028 IF (imin .GT. 1) THEN
1029 DO WHILE (phi(imin-1) .NE. zero)
1030 imin = imin - 1
1031 IF (imin .LE. 1) EXIT
1032 END DO
1033 END IF
1034*
1035* Repeat main iteration loop
1036*
1037 END DO
1038*
1039* Postprocessing: order THETA from least to greatest
1040*
1041 DO i = 1, q
1042*
1043 mini = i
1044 thetamin = theta(i)
1045 DO j = i+1, q
1046 IF( theta(j) .LT. thetamin ) THEN
1047 mini = j
1048 thetamin = theta(j)
1049 END IF
1050 END DO
1051*
1052 IF( mini .NE. i ) THEN
1053 theta(mini) = theta(i)
1054 theta(i) = thetamin
1055 IF( colmajor ) THEN
1056 IF( wantu1 )
1057 $ CALL cswap( p, u1(1,i), 1, u1(1,mini), 1 )
1058 IF( wantu2 )
1059 $ CALL cswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1060 IF( wantv1t )
1061 $ CALL cswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1062 IF( wantv2t )
1063 $ CALL cswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1064 $ ldv2t )
1065 ELSE
1066 IF( wantu1 )
1067 $ CALL cswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1068 IF( wantu2 )
1069 $ CALL cswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1070 IF( wantv1t )
1071 $ CALL cswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1072 IF( wantv2t )
1073 $ CALL cswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
1074 END IF
1075 END IF
1076*
1077 END DO
1078*
1079 RETURN
1080*
1081* End of CBBCSD
1082*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slartgp
subroutine slartgp(f, g, cs, sn, r)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition slartgp.f:95
slartgs
subroutine slartgs(x, y, sigma, cs, sn)
SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
Definition slartgs.f:90
slas2
subroutine slas2(f, g, h, ssmin, ssmax)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition slas2.f:105
clasr
subroutine clasr(side, pivot, direct, m, n, c, s, a, lda)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition clasr.f:200
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
cscal
subroutine cscal(n, ca, cx, incx)
CSCAL
Definition cscal.f:78
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: