LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ sbbcsd()

subroutine sbbcsd ( character  JOBU1,
character  JOBU2,
character  JOBV1T,
character  JOBV2T,
character  TRANS,
integer  M,
integer  P,
integer  Q,
real, dimension( * )  THETA,
real, dimension( * )  PHI,
real, dimension( ldu1, * )  U1,
integer  LDU1,
real, dimension( ldu2, * )  U2,
integer  LDU2,
real, dimension( ldv1t, * )  V1T,
integer  LDV1T,
real, 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( * )  WORK,
integer  LWORK,
integer  INFO 
)

SBBCSD

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

Purpose:
 SBBCSD 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 SORCSD 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 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (Q)
          When SBBCSD converges, B11D contains the cosines of THETA(1),
          ..., THETA(Q). If SBBCSD 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 SBBCSD converges, B11E contains zeros. If SBBCSD fails
          to converge, then B11E contains the superdiagonal of the
          partially reduced top-left block.
[out]B12D
          B12D is REAL array, dimension (Q)
          When SBBCSD converges, B12D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD 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 SBBCSD converges, B12E contains zeros. If SBBCSD fails
          to converge, then B12E contains the subdiagonal of the
          partially reduced top-right block.
[out]B21D
          B21D is REAL array, dimension (Q)
          When SBBCSD converges, B21D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD 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 SBBCSD converges, B21E contains zeros. If SBBCSD fails
          to converge, then B21E contains the subdiagonal of the
          partially reduced bottom-left block.
[out]B22D
          B22D is REAL array, dimension (Q)
          When SBBCSD converges, B22D contains the negative sines of
          THETA(1), ..., THETA(Q). If SBBCSD 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 SBBCSD converges, B22E contains zeros. If SBBCSD fails
          to converge, then B22E contains the subdiagonal of the
          partially reduced bottom-right block.
[out]WORK
          WORK is REAL 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 SBBCSD 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 sbbcsd.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  REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
343  $ B21D( * ), B21E( * ), B22D( * ), B22E( * ),
344  $ PHI( * ), THETA( * ), WORK( * )
345  REAL 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  REAL NEGONE
358  parameter( negone = -1.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  $ LWORKMIN, LWORKOPT, 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 slasr, sscal, sswap, 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 = 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( 'SBBCSD', -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  $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
589  ELSE
590  CALL slartgs( 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 slartgp( b11bulge, b11d(imin), work(iu1sn+imin-1),
618  $ work(iu1cs+imin-1), r )
619  ELSE IF( mu .LE. nu ) THEN
620  CALL slartgs( b11e( imin ), b11d( imin + 1 ), mu,
621  $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
622  ELSE
623  CALL slartgs( 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 slartgp( b21bulge, b21d(imin), work(iu2sn+imin-1),
628  $ work(iu2cs+imin-1), r )
629  ELSE IF( nu .LT. mu ) THEN
630  CALL slartgs( b21e( imin ), b21d( imin + 1 ), nu,
631  $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
632  ELSE
633  CALL slartgs( 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 slartgp( x2, x1, work(iv1tsn+i-1), work(iv1tcs+i-1),
701  $ r )
702  ELSE IF( .NOT. restart11 .AND. restart21 ) THEN
703  CALL slartgp( b11bulge, b11e(i-1), work(iv1tsn+i-1),
704  $ work(iv1tcs+i-1), r )
705  ELSE IF( restart11 .AND. .NOT. restart21 ) THEN
706  CALL slartgp( b21bulge, b21e(i-1), work(iv1tsn+i-1),
707  $ work(iv1tcs+i-1), r )
708  ELSE IF( mu .LE. nu ) THEN
709  CALL slartgs( b11d(i), b11e(i), mu, work(iv1tcs+i-1),
710  $ work(iv1tsn+i-1) )
711  ELSE
712  CALL slartgs( 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 slartgp( y2, y1, work(iv2tsn+i-1-1),
719  $ work(iv2tcs+i-1-1), r )
720  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
721  CALL slartgp( 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 slartgp( 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 slartgs( b12e(i-1), b12d(i), nu, work(iv2tcs+i-1-1),
728  $ work(iv2tsn+i-1-1) )
729  ELSE
730  CALL slartgs( 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 slartgp( x2, x1, work(iu1sn+i-1), work(iu1cs+i-1),
784  $ r )
785  ELSE IF( .NOT. restart11 .AND. restart12 ) THEN
786  CALL slartgp( b11bulge, b11d(i), work(iu1sn+i-1),
787  $ work(iu1cs+i-1), r )
788  ELSE IF( restart11 .AND. .NOT. restart12 ) THEN
789  CALL slartgp( b12bulge, b12e(i-1), work(iu1sn+i-1),
790  $ work(iu1cs+i-1), r )
791  ELSE IF( mu .LE. nu ) THEN
792  CALL slartgs( b11e(i), b11d(i+1), mu, work(iu1cs+i-1),
793  $ work(iu1sn+i-1) )
794  ELSE
795  CALL slartgs( 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 slartgp( y2, y1, work(iu2sn+i-1), work(iu2cs+i-1),
800  $ r )
801  ELSE IF( .NOT. restart21 .AND. restart22 ) THEN
802  CALL slartgp( b21bulge, b21d(i), work(iu2sn+i-1),
803  $ work(iu2cs+i-1), r )
804  ELSE IF( restart21 .AND. .NOT. restart22 ) THEN
805  CALL slartgp( b22bulge, b22e(i-1), work(iu2sn+i-1),
806  $ work(iu2cs+i-1), r )
807  ELSE IF( nu .LT. mu ) THEN
808  CALL slartgs( b21e(i), b21e(i+1), nu, work(iu2cs+i-1),
809  $ work(iu2sn+i-1) )
810  ELSE
811  CALL slartgs( 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 slartgp( y2, y1, work(iv2tsn+imax-1-1),
863  $ work(iv2tcs+imax-1-1), r )
864  ELSE IF( .NOT. restart12 .AND. restart22 ) THEN
865  CALL slartgp( 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 slartgp( 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 slartgs( b12e(imax-1), b12d(imax), nu,
872  $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
873  ELSE
874  CALL slartgs( 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 slasr( '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 slasr( '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 slasr( '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 slasr( '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 slasr( '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 slasr( '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 slasr( '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 slasr( '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 sscal( q, negone, v1t(imax,1), ldv1t )
944  ELSE
945  CALL sscal( 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 sscal( p, negone, u1(1,imax), 1 )
967  ELSE
968  CALL sscal( 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 sscal( m-p, negone, u2(1,imax), 1 )
977  ELSE
978  CALL sscal( 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 sscal( m-q, negone, v2t(imax,1), ldv2t )
989  ELSE
990  CALL sscal( 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 sswap( p, u1(1,i), 1, u1(1,mini), 1 )
1052  IF( wantu2 )
1053  $ CALL sswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1054  IF( wantv1t )
1055  $ CALL sswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1056  IF( wantv2t )
1057  $ CALL sswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1058  $ ldv2t )
1059  ELSE
1060  IF( wantu1 )
1061  $ CALL sswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1062  IF( wantu2 )
1063  $ CALL sswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1064  IF( wantv1t )
1065  $ CALL sswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1066  IF( wantv2t )
1067  $ CALL sswap( 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 SBBCSD
1076 *
slasr
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
Definition: slasr.f:199
slas2
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:107
slartgp
subroutine slartgp(F, G, CS, SN, R)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: slartgp.f:95
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
sswap
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
sscal
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: