LAPACK 3.12.1
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 326 of file cbbcsd.f.

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