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

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.

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 *
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
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
Definition: slas2.f:107
subroutine slartgp(F, G, CS, SN, R)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
Definition: slartgp.f:95
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
Definition: sswap.f:82
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
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: