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
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 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
DLASR applies a sequence of plane rotations to a general rectangular matrix.
logical function lsame(CA, CB)
LSAME
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 dlartgs(X, Y, SIGMA, CS, SN)
DLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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 dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlas2(F, G, H, SSMIN, SSMAX)
DLAS2 computes singular values of a 2-by-2 triangular matrix.