330 SUBROUTINE sbbcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q,
331 $ theta, phi, u1, ldu1, u2, ldu2, v1t, ldv1t,
332 $ v2t, ldv2t, b11d, b11e, b12d, b12e, b21d, b21e,
333 $ b22d, b22e, work, lwork, info )
341 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS
342 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LWORK, M, P, Q
345 REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
346 $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
347 $ phi( * ), theta( * ), work( * )
348 REAL U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
356 parameter ( maxitr = 6 )
357 REAL HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
358 parameter ( hundred = 100.0e0, meighth = -0.125e0,
359 $ one = 1.0e0, piover2 = 1.57079632679489662e0,
360 $ ten = 10.0e0, zero = 0.0e0 )
362 parameter ( negone = -1.0e0 )
365 LOGICAL COLMAJOR, LQUERY, RESTART11, RESTART12,
366 $ restart21, restart22, wantu1, wantu2, wantv1t,
368 INTEGER I, IMIN, IMAX, ITER, IU1CS, IU1SN, IU2CS,
369 $ iu2sn, iv1tcs, iv1tsn, iv2tcs, iv2tsn, j,
370 $ lworkmin, lworkopt, maxit, mini
371 REAL B11BULGE, B12BULGE, B21BULGE, B22BULGE, DUMMY,
372 $ eps, mu, nu, r, sigma11, sigma21,
373 $ temp, thetamax, thetamin, thresh, tol, tolmul,
374 $ unfl, x1, x2, y1, y2
383 EXTERNAL lsame, slamch
386 INTRINSIC abs, atan2, cos, max, min, sin, sqrt
393 lquery = lwork .EQ. -1
394 wantu1 = lsame( jobu1,
'Y' )
395 wantu2 = lsame( jobu2,
'Y' )
396 wantv1t = lsame( jobv1t,
'Y' )
397 wantv2t = lsame( jobv2t,
'Y' )
398 colmajor = .NOT. lsame( trans,
'T' )
402 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
404 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
406 ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q )
THEN
408 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
410 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
412 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
414 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
420 IF( info .EQ. 0 .AND. q .EQ. 0 )
THEN
428 IF( info .EQ. 0 )
THEN
437 lworkopt = iv2tsn + q - 1
440 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
445 IF( info .NE. 0 )
THEN
446 CALL xerbla(
'SBBCSD', -info )
448 ELSE IF( lquery )
THEN
454 eps = slamch(
'Epsilon' )
455 unfl = slamch(
'Safe minimum' )
456 tolmul = max( ten, min( hundred, eps**meighth ) )
458 thresh = max( tol, maxitr*q*q*unfl )
463 IF( theta(i) .LT. thresh )
THEN
465 ELSE IF( theta(i) .GT. piover2-thresh )
THEN
470 IF( phi(i) .LT. thresh )
THEN
472 ELSE IF( phi(i) .GT. piover2-thresh )
THEN
480 DO WHILE( imax .GT. 1 )
481 IF( phi(imax-1) .NE. zero )
THEN
487 IF ( imin .GT. 1 )
THEN
488 DO WHILE( phi(imin-1) .NE. zero )
490 IF ( imin .LE. 1 )
EXIT
501 DO WHILE( imax .GT. 1 )
505 b11d(imin) = cos( theta(imin) )
506 b21d(imin) = -sin( theta(imin) )
507 DO i = imin, imax - 1
508 b11e(i) = -sin( theta(i) ) * sin( phi(i) )
509 b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
510 b12d(i) = sin( theta(i) ) * cos( phi(i) )
511 b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
512 b21e(i) = -cos( theta(i) ) * sin( phi(i) )
513 b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
514 b22d(i) = cos( theta(i) ) * cos( phi(i) )
515 b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
517 b12d(imax) = sin( theta(imax) )
518 b22d(imax) = cos( theta(imax) )
522 IF( iter .GT. maxit )
THEN
525 IF( phi(i) .NE. zero )
531 iter = iter + imax - imin
535 thetamax = theta(imin)
536 thetamin = theta(imin)
538 IF( theta(i) > thetamax )
539 $ thetamax = theta(i)
540 IF( theta(i) < thetamin )
541 $ thetamin = theta(i)
544 IF( thetamax .GT. piover2 - thresh )
THEN
552 ELSE IF( thetamin .LT. thresh )
THEN
564 CALL slas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
566 CALL slas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
569 IF( sigma11 .LE. sigma21 )
THEN
571 nu = sqrt( one - mu**2 )
572 IF( mu .LT. thresh )
THEN
578 mu = sqrt( 1.0 - nu**2 )
579 IF( nu .LT. thresh )
THEN
588 IF( mu .LE. nu )
THEN
589 CALL slartgs( b11d(imin), b11e(imin), mu,
590 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
592 CALL slartgs( b21d(imin), b21e(imin), nu,
593 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
596 temp = work(iv1tcs+imin-1)*b11d(imin) +
597 $ work(iv1tsn+imin-1)*b11e(imin)
598 b11e(imin) = work(iv1tcs+imin-1)*b11e(imin) -
599 $ work(iv1tsn+imin-1)*b11d(imin)
601 b11bulge = work(iv1tsn+imin-1)*b11d(imin+1)
602 b11d(imin+1) = work(iv1tcs+imin-1)*b11d(imin+1)
603 temp = work(iv1tcs+imin-1)*b21d(imin) +
604 $ work(iv1tsn+imin-1)*b21e(imin)
605 b21e(imin) = work(iv1tcs+imin-1)*b21e(imin) -
606 $ work(iv1tsn+imin-1)*b21d(imin)
608 b21bulge = work(iv1tsn+imin-1)*b21d(imin+1)
609 b21d(imin+1) = work(iv1tcs+imin-1)*b21d(imin+1)
613 theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
614 $ sqrt( b11d(imin)**2+b11bulge**2 ) )
618 IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 )
THEN
619 CALL slartgp( b11bulge, b11d(imin), work(iu1sn+imin-1),
620 $ work(iu1cs+imin-1), r )
621 ELSE IF( mu .LE. nu )
THEN
622 CALL slartgs( b11e( imin ), b11d( imin + 1 ), mu,
623 $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
625 CALL slartgs( b12d( imin ), b12e( imin ), nu,
626 $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
628 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 )
THEN
629 CALL slartgp( b21bulge, b21d(imin), work(iu2sn+imin-1),
630 $ work(iu2cs+imin-1), r )
631 ELSE IF( nu .LT. mu )
THEN
632 CALL slartgs( b21e( imin ), b21d( imin + 1 ), nu,
633 $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
635 CALL slartgs( b22d(imin), b22e(imin), mu,
636 $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
638 work(iu2cs+imin-1) = -work(iu2cs+imin-1)
639 work(iu2sn+imin-1) = -work(iu2sn+imin-1)
641 temp = work(iu1cs+imin-1)*b11e(imin) +
642 $ work(iu1sn+imin-1)*b11d(imin+1)
643 b11d(imin+1) = work(iu1cs+imin-1)*b11d(imin+1) -
644 $ work(iu1sn+imin-1)*b11e(imin)
646 IF( imax .GT. imin+1 )
THEN
647 b11bulge = work(iu1sn+imin-1)*b11e(imin+1)
648 b11e(imin+1) = work(iu1cs+imin-1)*b11e(imin+1)
650 temp = work(iu1cs+imin-1)*b12d(imin) +
651 $ work(iu1sn+imin-1)*b12e(imin)
652 b12e(imin) = work(iu1cs+imin-1)*b12e(imin) -
653 $ work(iu1sn+imin-1)*b12d(imin)
655 b12bulge = work(iu1sn+imin-1)*b12d(imin+1)
656 b12d(imin+1) = work(iu1cs+imin-1)*b12d(imin+1)
657 temp = work(iu2cs+imin-1)*b21e(imin) +
658 $ work(iu2sn+imin-1)*b21d(imin+1)
659 b21d(imin+1) = work(iu2cs+imin-1)*b21d(imin+1) -
660 $ work(iu2sn+imin-1)*b21e(imin)
662 IF( imax .GT. imin+1 )
THEN
663 b21bulge = work(iu2sn+imin-1)*b21e(imin+1)
664 b21e(imin+1) = work(iu2cs+imin-1)*b21e(imin+1)
666 temp = work(iu2cs+imin-1)*b22d(imin) +
667 $ work(iu2sn+imin-1)*b22e(imin)
668 b22e(imin) = work(iu2cs+imin-1)*b22e(imin) -
669 $ work(iu2sn+imin-1)*b22d(imin)
671 b22bulge = work(iu2sn+imin-1)*b22d(imin+1)
672 b22d(imin+1) = work(iu2cs+imin-1)*b22d(imin+1)
678 DO i = imin+1, imax-1
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
687 phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
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
701 IF( .NOT. restart11 .AND. .NOT. restart21 )
THEN
702 CALL slartgp( x2, x1, work(iv1tsn+i-1), work(iv1tcs+i-1),
704 ELSE IF( .NOT. restart11 .AND. restart21 )
THEN
705 CALL slartgp( b11bulge, b11e(i-1), work(iv1tsn+i-1),
706 $ work(iv1tcs+i-1), r )
707 ELSE IF( restart11 .AND. .NOT. restart21 )
THEN
708 CALL slartgp( b21bulge, b21e(i-1), work(iv1tsn+i-1),
709 $ work(iv1tcs+i-1), r )
710 ELSE IF( mu .LE. nu )
THEN
711 CALL slartgs( b11d(i), b11e(i), mu, work(iv1tcs+i-1),
714 CALL slartgs( b21d(i), b21e(i), nu, work(iv1tcs+i-1),
717 work(iv1tcs+i-1) = -work(iv1tcs+i-1)
718 work(iv1tsn+i-1) = -work(iv1tsn+i-1)
719 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
720 CALL slartgp( y2, y1, work(iv2tsn+i-1-1),
721 $ work(iv2tcs+i-1-1), r )
722 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
723 CALL slartgp( b12bulge, b12d(i-1), work(iv2tsn+i-1-1),
724 $ work(iv2tcs+i-1-1), r )
725 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
726 CALL slartgp( b22bulge, b22d(i-1), work(iv2tsn+i-1-1),
727 $ work(iv2tcs+i-1-1), r )
728 ELSE IF( nu .LT. mu )
THEN
729 CALL slartgs( b12e(i-1), b12d(i), nu, work(iv2tcs+i-1-1),
730 $ work(iv2tsn+i-1-1) )
732 CALL slartgs( b22e(i-1), b22d(i), mu, work(iv2tcs+i-1-1),
733 $ work(iv2tsn+i-1-1) )
736 temp = work(iv1tcs+i-1)*b11d(i) + work(iv1tsn+i-1)*b11e(i)
737 b11e(i) = work(iv1tcs+i-1)*b11e(i) -
738 $ work(iv1tsn+i-1)*b11d(i)
740 b11bulge = work(iv1tsn+i-1)*b11d(i+1)
741 b11d(i+1) = work(iv1tcs+i-1)*b11d(i+1)
742 temp = work(iv1tcs+i-1)*b21d(i) + work(iv1tsn+i-1)*b21e(i)
743 b21e(i) = work(iv1tcs+i-1)*b21e(i) -
744 $ work(iv1tsn+i-1)*b21d(i)
746 b21bulge = work(iv1tsn+i-1)*b21d(i+1)
747 b21d(i+1) = work(iv1tcs+i-1)*b21d(i+1)
748 temp = work(iv2tcs+i-1-1)*b12e(i-1) +
749 $ work(iv2tsn+i-1-1)*b12d(i)
750 b12d(i) = work(iv2tcs+i-1-1)*b12d(i) -
751 $ work(iv2tsn+i-1-1)*b12e(i-1)
753 b12bulge = work(iv2tsn+i-1-1)*b12e(i)
754 b12e(i) = work(iv2tcs+i-1-1)*b12e(i)
755 temp = work(iv2tcs+i-1-1)*b22e(i-1) +
756 $ work(iv2tsn+i-1-1)*b22d(i)
757 b22d(i) = work(iv2tcs+i-1-1)*b22d(i) -
758 $ work(iv2tsn+i-1-1)*b22e(i-1)
760 b22bulge = work(iv2tsn+i-1-1)*b22e(i)
761 b22e(i) = work(iv2tcs+i-1-1)*b22e(i)
765 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
766 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
767 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
768 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
770 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
775 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
776 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
777 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
778 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
784 IF( .NOT. restart11 .AND. .NOT. restart12 )
THEN
785 CALL slartgp( x2, x1, work(iu1sn+i-1), work(iu1cs+i-1),
787 ELSE IF( .NOT. restart11 .AND. restart12 )
THEN
788 CALL slartgp( b11bulge, b11d(i), work(iu1sn+i-1),
789 $ work(iu1cs+i-1), r )
790 ELSE IF( restart11 .AND. .NOT. restart12 )
THEN
791 CALL slartgp( b12bulge, b12e(i-1), work(iu1sn+i-1),
792 $ work(iu1cs+i-1), r )
793 ELSE IF( mu .LE. nu )
THEN
794 CALL slartgs( b11e(i), b11d(i+1), mu, work(iu1cs+i-1),
797 CALL slartgs( b12d(i), b12e(i), nu, work(iu1cs+i-1),
800 IF( .NOT. restart21 .AND. .NOT. restart22 )
THEN
801 CALL slartgp( y2, y1, work(iu2sn+i-1), work(iu2cs+i-1),
803 ELSE IF( .NOT. restart21 .AND. restart22 )
THEN
804 CALL slartgp( b21bulge, b21d(i), work(iu2sn+i-1),
805 $ work(iu2cs+i-1), r )
806 ELSE IF( restart21 .AND. .NOT. restart22 )
THEN
807 CALL slartgp( b22bulge, b22e(i-1), work(iu2sn+i-1),
808 $ work(iu2cs+i-1), r )
809 ELSE IF( nu .LT. mu )
THEN
810 CALL slartgs( b21e(i), b21e(i+1), nu, work(iu2cs+i-1),
813 CALL slartgs( b22d(i), b22e(i), mu, work(iu2cs+i-1),
816 work(iu2cs+i-1) = -work(iu2cs+i-1)
817 work(iu2sn+i-1) = -work(iu2sn+i-1)
819 temp = work(iu1cs+i-1)*b11e(i) + work(iu1sn+i-1)*b11d(i+1)
820 b11d(i+1) = work(iu1cs+i-1)*b11d(i+1) -
821 $ work(iu1sn+i-1)*b11e(i)
823 IF( i .LT. imax - 1 )
THEN
824 b11bulge = work(iu1sn+i-1)*b11e(i+1)
825 b11e(i+1) = work(iu1cs+i-1)*b11e(i+1)
827 temp = work(iu2cs+i-1)*b21e(i) + work(iu2sn+i-1)*b21d(i+1)
828 b21d(i+1) = work(iu2cs+i-1)*b21d(i+1) -
829 $ work(iu2sn+i-1)*b21e(i)
831 IF( i .LT. imax - 1 )
THEN
832 b21bulge = work(iu2sn+i-1)*b21e(i+1)
833 b21e(i+1) = work(iu2cs+i-1)*b21e(i+1)
835 temp = work(iu1cs+i-1)*b12d(i) + work(iu1sn+i-1)*b12e(i)
836 b12e(i) = work(iu1cs+i-1)*b12e(i) - work(iu1sn+i-1)*b12d(i)
838 b12bulge = work(iu1sn+i-1)*b12d(i+1)
839 b12d(i+1) = work(iu1cs+i-1)*b12d(i+1)
840 temp = work(iu2cs+i-1)*b22d(i) + work(iu2sn+i-1)*b22e(i)
841 b22e(i) = work(iu2cs+i-1)*b22e(i) - work(iu2sn+i-1)*b22d(i)
843 b22bulge = work(iu2sn+i-1)*b22d(i+1)
844 b22d(i+1) = work(iu2cs+i-1)*b22d(i+1)
850 x1 = sin(theta(imax-1))*b11e(imax-1) +
851 $ cos(theta(imax-1))*b21e(imax-1)
852 y1 = sin(theta(imax-1))*b12d(imax-1) +
853 $ cos(theta(imax-1))*b22d(imax-1)
854 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
856 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
860 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
861 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
863 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
864 CALL slartgp( y2, y1, work(iv2tsn+imax-1-1),
865 $ work(iv2tcs+imax-1-1), r )
866 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
867 CALL slartgp( b12bulge, b12d(imax-1), work(iv2tsn+imax-1-1),
868 $ work(iv2tcs+imax-1-1), r )
869 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
870 CALL slartgp( b22bulge, b22d(imax-1), work(iv2tsn+imax-1-1),
871 $ work(iv2tcs+imax-1-1), r )
872 ELSE IF( nu .LT. mu )
THEN
873 CALL slartgs( b12e(imax-1), b12d(imax), nu,
874 $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
876 CALL slartgs( b22e(imax-1), b22d(imax), mu,
877 $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
880 temp = work(iv2tcs+imax-1-1)*b12e(imax-1) +
881 $ work(iv2tsn+imax-1-1)*b12d(imax)
882 b12d(imax) = work(iv2tcs+imax-1-1)*b12d(imax) -
883 $ work(iv2tsn+imax-1-1)*b12e(imax-1)
885 temp = work(iv2tcs+imax-1-1)*b22e(imax-1) +
886 $ work(iv2tsn+imax-1-1)*b22d(imax)
887 b22d(imax) = work(iv2tcs+imax-1-1)*b22d(imax) -
888 $ work(iv2tsn+imax-1-1)*b22e(imax-1)
895 CALL slasr(
'R',
'V',
'F', p, imax-imin+1,
896 $ work(iu1cs+imin-1), work(iu1sn+imin-1),
899 CALL slasr(
'L',
'V',
'F', imax-imin+1, p,
900 $ work(iu1cs+imin-1), work(iu1sn+imin-1),
906 CALL slasr(
'R',
'V',
'F', m-p, imax-imin+1,
907 $ work(iu2cs+imin-1), work(iu2sn+imin-1),
910 CALL slasr(
'L',
'V',
'F', imax-imin+1, m-p,
911 $ work(iu2cs+imin-1), work(iu2sn+imin-1),
917 CALL slasr(
'L',
'V',
'F', imax-imin+1, q,
918 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
919 $ v1t(imin,1), ldv1t )
921 CALL slasr(
'R',
'V',
'F', q, imax-imin+1,
922 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
923 $ v1t(1,imin), ldv1t )
928 CALL slasr(
'L',
'V',
'F', imax-imin+1, m-q,
929 $ work(iv2tcs+imin-1), work(iv2tsn+imin-1),
930 $ v2t(imin,1), ldv2t )
932 CALL slasr(
'R',
'V',
'F', m-q, imax-imin+1,
933 $ work(iv2tcs+imin-1), work(iv2tsn+imin-1),
934 $ v2t(1,imin), ldv2t )
940 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 )
THEN
941 b11d(imax) = -b11d(imax)
942 b21d(imax) = -b21d(imax)
945 CALL sscal( q, negone, v1t(imax,1), ldv1t )
947 CALL sscal( q, negone, v1t(1,imax), 1 )
954 x1 = cos(phi(imax-1))*b11d(imax) +
955 $ sin(phi(imax-1))*b12e(imax-1)
956 y1 = cos(phi(imax-1))*b21d(imax) +
957 $ sin(phi(imax-1))*b22e(imax-1)
959 theta(imax) = atan2( abs(y1), abs(x1) )
964 IF( b11d(imax)+b12e(imax-1) .LT. 0 )
THEN
965 b12d(imax) = -b12d(imax)
968 CALL sscal( p, negone, u1(1,imax), 1 )
970 CALL sscal( p, negone, u1(imax,1), ldu1 )
974 IF( b21d(imax)+b22e(imax-1) .GT. 0 )
THEN
975 b22d(imax) = -b22d(imax)
978 CALL sscal( m-p, negone, u2(1,imax), 1 )
980 CALL sscal( m-p, negone, u2(imax,1), ldu2 )
987 IF( b12d(imax)+b22d(imax) .LT. 0 )
THEN
990 CALL sscal( m-q, negone, v2t(imax,1), ldv2t )
992 CALL sscal( m-q, negone, v2t(1,imax), 1 )
1000 IF( theta(i) .LT. thresh )
THEN
1002 ELSE IF( theta(i) .GT. piover2-thresh )
THEN
1007 IF( phi(i) .LT. thresh )
THEN
1009 ELSE IF( phi(i) .GT. piover2-thresh )
THEN
1016 IF (imax .GT. 1)
THEN
1017 DO WHILE( phi(imax-1) .EQ. zero )
1019 IF (imax .LE. 1)
EXIT
1022 IF( imin .GT. imax - 1 )
1024 IF (imin .GT. 1)
THEN
1025 DO WHILE (phi(imin-1) .NE. zero)
1027 IF (imin .LE. 1)
EXIT
1042 IF( theta(j) .LT. thetamin )
THEN
1048 IF( mini .NE. i )
THEN
1049 theta(mini) = theta(i)
1053 $
CALL sswap( p, u1(1,i), 1, u1(1,mini), 1 )
1055 $
CALL sswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1057 $
CALL sswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1059 $
CALL sswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1063 $
CALL sswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1065 $
CALL sswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1067 $
CALL sswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1069 $
CALL sswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
subroutine sbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, WORK, LWORK, INFO)
SBBCSD
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
subroutine slasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
SLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine slartgp(F, G, CS, SN, R)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine slartgs(X, Y, SIGMA, CS, SN)
SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...