330 SUBROUTINE dbbcsd( 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 DOUBLE PRECISION B11D( * ), B11E( * ), B12D( * ), B12E( * ),
346 $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
347 $ phi( * ), theta( * ), work( * )
348 DOUBLE PRECISION U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
356 parameter ( maxitr = 6 )
357 DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
358 parameter ( hundred = 100.0d0, meighth = -0.125d0,
359 $ one = 1.0d0, piover2 = 1.57079632679489662d0,
360 $ ten = 10.0d0, zero = 0.0d0 )
361 DOUBLE PRECISION NEGONE
362 parameter ( negone = -1.0d0 )
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 DOUBLE PRECISION 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
381 DOUBLE PRECISION DLAMCH
383 EXTERNAL lsame, dlamch
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(
'DBBCSD', -info )
448 ELSE IF( lquery )
THEN
454 eps = dlamch(
'Epsilon' )
455 unfl = dlamch(
'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 dlas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
566 CALL dlas2( 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 dlartgs( b11d(imin), b11e(imin), mu,
590 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1) )
592 CALL dlartgs( 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 dlartgp( b11bulge, b11d(imin), work(iu1sn+imin-1),
620 $ work(iu1cs+imin-1), r )
621 ELSE IF( mu .LE. nu )
THEN
622 CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
623 $ work(iu1cs+imin-1), work(iu1sn+imin-1) )
625 CALL dlartgs( 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 dlartgp( b21bulge, b21d(imin), work(iu2sn+imin-1),
630 $ work(iu2cs+imin-1), r )
631 ELSE IF( nu .LT. mu )
THEN
632 CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
633 $ work(iu2cs+imin-1), work(iu2sn+imin-1) )
635 CALL dlartgs( 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 dlartgp( x2, x1, work(iv1tsn+i-1), work(iv1tcs+i-1),
704 ELSE IF( .NOT. restart11 .AND. restart21 )
THEN
705 CALL dlartgp( b11bulge, b11e(i-1), work(iv1tsn+i-1),
706 $ work(iv1tcs+i-1), r )
707 ELSE IF( restart11 .AND. .NOT. restart21 )
THEN
708 CALL dlartgp( b21bulge, b21e(i-1), work(iv1tsn+i-1),
709 $ work(iv1tcs+i-1), r )
710 ELSE IF( mu .LE. nu )
THEN
711 CALL dlartgs( b11d(i), b11e(i), mu, work(iv1tcs+i-1),
714 CALL dlartgs( 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 dlartgp( y2, y1, work(iv2tsn+i-1-1),
721 $ work(iv2tcs+i-1-1), r )
722 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
723 CALL dlartgp( 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 dlartgp( 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 dlartgs( b12e(i-1), b12d(i), nu, work(iv2tcs+i-1-1),
730 $ work(iv2tsn+i-1-1) )
732 CALL dlartgs( 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 dlartgp( x2, x1, work(iu1sn+i-1), work(iu1cs+i-1),
787 ELSE IF( .NOT. restart11 .AND. restart12 )
THEN
788 CALL dlartgp( b11bulge, b11d(i), work(iu1sn+i-1),
789 $ work(iu1cs+i-1), r )
790 ELSE IF( restart11 .AND. .NOT. restart12 )
THEN
791 CALL dlartgp( b12bulge, b12e(i-1), work(iu1sn+i-1),
792 $ work(iu1cs+i-1), r )
793 ELSE IF( mu .LE. nu )
THEN
794 CALL dlartgs( b11e(i), b11d(i+1), mu, work(iu1cs+i-1),
797 CALL dlartgs( b12d(i), b12e(i), nu, work(iu1cs+i-1),
800 IF( .NOT. restart21 .AND. .NOT. restart22 )
THEN
801 CALL dlartgp( y2, y1, work(iu2sn+i-1), work(iu2cs+i-1),
803 ELSE IF( .NOT. restart21 .AND. restart22 )
THEN
804 CALL dlartgp( b21bulge, b21d(i), work(iu2sn+i-1),
805 $ work(iu2cs+i-1), r )
806 ELSE IF( restart21 .AND. .NOT. restart22 )
THEN
807 CALL dlartgp( b22bulge, b22e(i-1), work(iu2sn+i-1),
808 $ work(iu2cs+i-1), r )
809 ELSE IF( nu .LT. mu )
THEN
810 CALL dlartgs( b21e(i), b21e(i+1), nu, work(iu2cs+i-1),
813 CALL dlartgs( 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 dlartgp( y2, y1, work(iv2tsn+imax-1-1),
865 $ work(iv2tcs+imax-1-1), r )
866 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
867 CALL dlartgp( 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 dlartgp( 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 dlartgs( b12e(imax-1), b12d(imax), nu,
874 $ work(iv2tcs+imax-1-1), work(iv2tsn+imax-1-1) )
876 CALL dlartgs( 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 dlasr(
'R',
'V',
'F', p, imax-imin+1,
896 $ work(iu1cs+imin-1), work(iu1sn+imin-1),
899 CALL dlasr(
'L',
'V',
'F', imax-imin+1, p,
900 $ work(iu1cs+imin-1), work(iu1sn+imin-1),
906 CALL dlasr(
'R',
'V',
'F', m-p, imax-imin+1,
907 $ work(iu2cs+imin-1), work(iu2sn+imin-1),
910 CALL dlasr(
'L',
'V',
'F', imax-imin+1, m-p,
911 $ work(iu2cs+imin-1), work(iu2sn+imin-1),
917 CALL dlasr(
'L',
'V',
'F', imax-imin+1, q,
918 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
919 $ v1t(imin,1), ldv1t )
921 CALL dlasr(
'R',
'V',
'F', q, imax-imin+1,
922 $ work(iv1tcs+imin-1), work(iv1tsn+imin-1),
923 $ v1t(1,imin), ldv1t )
928 CALL dlasr(
'L',
'V',
'F', imax-imin+1, m-q,
929 $ work(iv2tcs+imin-1), work(iv2tsn+imin-1),
930 $ v2t(imin,1), ldv2t )
932 CALL dlasr(
'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 dscal( q, negone, v1t(imax,1), ldv1t )
947 CALL dscal( 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 dscal( p, negone, u1(1,imax), 1 )
970 CALL dscal( p, negone, u1(imax,1), ldu1 )
974 IF( b21d(imax)+b22e(imax-1) .GT. 0 )
THEN
975 b22d(imax) = -b22d(imax)
978 CALL dscal( m-p, negone, u2(1,imax), 1 )
980 CALL dscal( m-p, negone, u2(imax,1), ldu2 )
987 IF( b12d(imax)+b22d(imax) .LT. 0 )
THEN
990 CALL dscal( m-q, negone, v2t(imax,1), ldv2t )
992 CALL dscal( 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 dswap( p, u1(1,i), 1, u1(1,mini), 1 )
1055 $
CALL dswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1057 $
CALL dswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1059 $
CALL dswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1063 $
CALL dswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1065 $
CALL dswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1067 $
CALL dswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1069 $
CALL dswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine dlartgs(X, Y, SIGMA, CS, SN)
DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
subroutine dscal(N, DA, DX, INCX)
DSCAL
subroutine dlartgp(F, G, CS, SN, R)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
subroutine dbbcsd(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)
DBBCSD
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.