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
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
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 )
double precision function dlamch(CMACH)
DLAMCH
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.
logical function lsame(CA, CB)
LSAME
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.