341 CHARACTER jobu1, jobu2, jobv1t, jobv2t, trans
342 INTEGER info, ldu1, ldu2, ldv1t, ldv2t, lrwork, m, p, q
345 DOUBLE PRECISION b11d( * ), b11e( * ), b12d( * ), b12e( * ),
346 $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
347 $ phi( * ), theta( * ), rwork( * )
348 COMPLEX*16 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 COMPLEX*16 negonecomplex
362 parameter ( negonecomplex = (-1.0d0,0.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 $ lrworkmin, lrworkopt, 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
385 INTRINSIC abs, atan2, cos, max, min, sin, sqrt
392 lquery = lrwork .EQ. -1
393 wantu1 =
lsame( jobu1,
'Y' )
394 wantu2 =
lsame( jobu2,
'Y' )
395 wantv1t =
lsame( jobv1t,
'Y' )
396 wantv2t =
lsame( jobv2t,
'Y' )
397 colmajor = .NOT.
lsame( trans,
'T' )
401 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
403 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
405 ELSE IF( q .GT. p .OR. q .GT. m-p .OR. q .GT. m-q )
THEN
407 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
409 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
411 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
413 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
419 IF( info .EQ. 0 .AND. q .EQ. 0 )
THEN
427 IF( info .EQ. 0 )
THEN
436 lrworkopt = iv2tsn + q - 1
437 lrworkmin = lrworkopt
439 IF( lrwork .LT. lrworkmin .AND. .NOT. lquery )
THEN
444 IF( info .NE. 0 )
THEN
445 CALL xerbla(
'ZBBCSD', -info )
447 ELSE IF( lquery )
THEN
454 unfl =
dlamch(
'Safe minimum' )
455 tolmul = max( ten, min( hundred, eps**meighth ) )
457 thresh = max( tol, maxitr*q*q*unfl )
462 IF( theta(i) .LT. thresh )
THEN
464 ELSE IF( theta(i) .GT. piover2-thresh )
THEN
469 IF( phi(i) .LT. thresh )
THEN
471 ELSE IF( phi(i) .GT. piover2-thresh )
THEN
479 DO WHILE( imax .GT. 1 )
480 IF( phi(imax-1) .NE. zero )
THEN
486 IF ( imin .GT. 1 )
THEN
487 DO WHILE( phi(imin-1) .NE. zero )
489 IF ( imin .LE. 1 )
EXIT
500 DO WHILE( imax .GT. 1 )
504 b11d(imin) = cos( theta(imin) )
505 b21d(imin) = -sin( theta(imin) )
506 DO i = imin, imax - 1
507 b11e(i) = -sin( theta(i) ) * sin( phi(i) )
508 b11d(i+1) = cos( theta(i+1) ) * cos( phi(i) )
509 b12d(i) = sin( theta(i) ) * cos( phi(i) )
510 b12e(i) = cos( theta(i+1) ) * sin( phi(i) )
511 b21e(i) = -cos( theta(i) ) * sin( phi(i) )
512 b21d(i+1) = -sin( theta(i+1) ) * cos( phi(i) )
513 b22d(i) = cos( theta(i) ) * cos( phi(i) )
514 b22e(i) = -sin( theta(i+1) ) * sin( phi(i) )
516 b12d(imax) = sin( theta(imax) )
517 b22d(imax) = cos( theta(imax) )
521 IF( iter .GT. maxit )
THEN
524 IF( phi(i) .NE. zero )
530 iter = iter + imax - imin
534 thetamax = theta(imin)
535 thetamin = theta(imin)
537 IF( theta(i) > thetamax )
538 $ thetamax = theta(i)
539 IF( theta(i) < thetamin )
540 $ thetamin = theta(i)
543 IF( thetamax .GT. piover2 - thresh )
THEN
551 ELSE IF( thetamin .LT. thresh )
THEN
563 CALL dlas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
565 CALL dlas2( b21d(imax-1), b21e(imax-1), b21d(imax), sigma21,
568 IF( sigma11 .LE. sigma21 )
THEN
570 nu = sqrt( one - mu**2 )
571 IF( mu .LT. thresh )
THEN
577 mu = sqrt( 1.0 - nu**2 )
578 IF( nu .LT. thresh )
THEN
587 IF( mu .LE. nu )
THEN
588 CALL dlartgs( b11d(imin), b11e(imin), mu,
589 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
591 CALL dlartgs( b21d(imin), b21e(imin), nu,
592 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
595 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
596 $ rwork(iv1tsn+imin-1)*b11e(imin)
597 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
598 $ rwork(iv1tsn+imin-1)*b11d(imin)
600 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
601 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
602 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
603 $ rwork(iv1tsn+imin-1)*b21e(imin)
604 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
605 $ rwork(iv1tsn+imin-1)*b21d(imin)
607 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
608 b21d(imin+1) = rwork(iv1tcs+imin-1)*b21d(imin+1)
612 theta( imin ) = atan2( sqrt( b21d(imin)**2+b21bulge**2 ),
613 $ sqrt( b11d(imin)**2+b11bulge**2 ) )
617 IF( b11d(imin)**2+b11bulge**2 .GT. thresh**2 )
THEN
618 CALL dlartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
619 $ rwork(iu1cs+imin-1), r )
620 ELSE IF( mu .LE. nu )
THEN
621 CALL dlartgs( b11e( imin ), b11d( imin + 1 ), mu,
622 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
624 CALL dlartgs( b12d( imin ), b12e( imin ), nu,
625 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
627 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 )
THEN
628 CALL dlartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
629 $ rwork(iu2cs+imin-1), r )
630 ELSE IF( nu .LT. mu )
THEN
631 CALL dlartgs( b21e( imin ), b21d( imin + 1 ), nu,
632 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
634 CALL dlartgs( b22d(imin), b22e(imin), mu,
635 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
637 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
638 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
640 temp = rwork(iu1cs+imin-1)*b11e(imin) +
641 $ rwork(iu1sn+imin-1)*b11d(imin+1)
642 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
643 $ rwork(iu1sn+imin-1)*b11e(imin)
645 IF( imax .GT. imin+1 )
THEN
646 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
647 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
649 temp = rwork(iu1cs+imin-1)*b12d(imin) +
650 $ rwork(iu1sn+imin-1)*b12e(imin)
651 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
652 $ rwork(iu1sn+imin-1)*b12d(imin)
654 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
655 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
656 temp = rwork(iu2cs+imin-1)*b21e(imin) +
657 $ rwork(iu2sn+imin-1)*b21d(imin+1)
658 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
659 $ rwork(iu2sn+imin-1)*b21e(imin)
661 IF( imax .GT. imin+1 )
THEN
662 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
663 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
665 temp = rwork(iu2cs+imin-1)*b22d(imin) +
666 $ rwork(iu2sn+imin-1)*b22e(imin)
667 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
668 $ rwork(iu2sn+imin-1)*b22d(imin)
670 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
671 b22d(imin+1) = rwork(iu2cs+imin-1)*b22d(imin+1)
677 DO i = imin+1, imax-1
681 x1 = sin(theta(i-1))*b11e(i-1) + cos(theta(i-1))*b21e(i-1)
682 x2 = sin(theta(i-1))*b11bulge + cos(theta(i-1))*b21bulge
683 y1 = sin(theta(i-1))*b12d(i-1) + cos(theta(i-1))*b22d(i-1)
684 y2 = sin(theta(i-1))*b12bulge + cos(theta(i-1))*b22bulge
686 phi(i-1) = atan2( sqrt(x1**2+x2**2), sqrt(y1**2+y2**2) )
691 restart11 = b11e(i-1)**2 + b11bulge**2 .LE. thresh**2
692 restart21 = b21e(i-1)**2 + b21bulge**2 .LE. thresh**2
693 restart12 = b12d(i-1)**2 + b12bulge**2 .LE. thresh**2
694 restart22 = b22d(i-1)**2 + b22bulge**2 .LE. thresh**2
700 IF( .NOT. restart11 .AND. .NOT. restart21 )
THEN
701 CALL dlartgp( x2, x1, rwork(iv1tsn+i-1),
702 $ rwork(iv1tcs+i-1), r )
703 ELSE IF( .NOT. restart11 .AND. restart21 )
THEN
704 CALL dlartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
705 $ rwork(iv1tcs+i-1), r )
706 ELSE IF( restart11 .AND. .NOT. restart21 )
THEN
707 CALL dlartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
708 $ rwork(iv1tcs+i-1), r )
709 ELSE IF( mu .LE. nu )
THEN
710 CALL dlartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
711 $ rwork(iv1tsn+i-1) )
713 CALL dlartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
714 $ rwork(iv1tsn+i-1) )
716 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
717 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
718 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
719 CALL dlartgp( y2, y1, rwork(iv2tsn+i-1-1),
720 $ rwork(iv2tcs+i-1-1), r )
721 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
722 CALL dlartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
723 $ rwork(iv2tcs+i-1-1), r )
724 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
725 CALL dlartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
726 $ rwork(iv2tcs+i-1-1), r )
727 ELSE IF( nu .LT. mu )
THEN
728 CALL dlartgs( b12e(i-1), b12d(i), nu,
729 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
731 CALL dlartgs( b22e(i-1), b22d(i), mu,
732 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
735 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
736 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
737 $ rwork(iv1tsn+i-1)*b11d(i)
739 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
740 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
741 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
742 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
743 $ rwork(iv1tsn+i-1)*b21d(i)
745 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
746 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
747 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
748 $ rwork(iv2tsn+i-1-1)*b12d(i)
749 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
750 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
752 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
753 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
754 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
755 $ rwork(iv2tsn+i-1-1)*b22d(i)
756 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
757 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
759 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
760 b22e(i) = rwork(iv2tcs+i-1-1)*b22e(i)
764 x1 = cos(phi(i-1))*b11d(i) + sin(phi(i-1))*b12e(i-1)
765 x2 = cos(phi(i-1))*b11bulge + sin(phi(i-1))*b12bulge
766 y1 = cos(phi(i-1))*b21d(i) + sin(phi(i-1))*b22e(i-1)
767 y2 = cos(phi(i-1))*b21bulge + sin(phi(i-1))*b22bulge
769 theta(i) = atan2( sqrt(y1**2+y2**2), sqrt(x1**2+x2**2) )
774 restart11 = b11d(i)**2 + b11bulge**2 .LE. thresh**2
775 restart12 = b12e(i-1)**2 + b12bulge**2 .LE. thresh**2
776 restart21 = b21d(i)**2 + b21bulge**2 .LE. thresh**2
777 restart22 = b22e(i-1)**2 + b22bulge**2 .LE. thresh**2
783 IF( .NOT. restart11 .AND. .NOT. restart12 )
THEN
784 CALL dlartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
786 ELSE IF( .NOT. restart11 .AND. restart12 )
THEN
787 CALL dlartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
788 $ rwork(iu1cs+i-1), r )
789 ELSE IF( restart11 .AND. .NOT. restart12 )
THEN
790 CALL dlartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
791 $ rwork(iu1cs+i-1), r )
792 ELSE IF( mu .LE. nu )
THEN
793 CALL dlartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
796 CALL dlartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
799 IF( .NOT. restart21 .AND. .NOT. restart22 )
THEN
800 CALL dlartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
802 ELSE IF( .NOT. restart21 .AND. restart22 )
THEN
803 CALL dlartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
804 $ rwork(iu2cs+i-1), r )
805 ELSE IF( restart21 .AND. .NOT. restart22 )
THEN
806 CALL dlartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
807 $ rwork(iu2cs+i-1), r )
808 ELSE IF( nu .LT. mu )
THEN
809 CALL dlartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
812 CALL dlartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
815 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
816 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
818 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
819 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
820 $ rwork(iu1sn+i-1)*b11e(i)
822 IF( i .LT. imax - 1 )
THEN
823 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
824 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
826 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
827 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
828 $ rwork(iu2sn+i-1)*b21e(i)
830 IF( i .LT. imax - 1 )
THEN
831 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
832 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
834 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
835 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
836 $ rwork(iu1sn+i-1)*b12d(i)
838 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
839 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
840 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
841 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
842 $ rwork(iu2sn+i-1)*b22d(i)
844 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
845 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
851 x1 = sin(theta(imax-1))*b11e(imax-1) +
852 $ cos(theta(imax-1))*b21e(imax-1)
853 y1 = sin(theta(imax-1))*b12d(imax-1) +
854 $ cos(theta(imax-1))*b22d(imax-1)
855 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
857 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
861 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
862 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
864 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
865 CALL dlartgp( y2, y1, rwork(iv2tsn+imax-1-1),
866 $ rwork(iv2tcs+imax-1-1), r )
867 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
868 CALL dlartgp( b12bulge, b12d(imax-1),
869 $ rwork(iv2tsn+imax-1-1),
870 $ rwork(iv2tcs+imax-1-1), r )
871 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
872 CALL dlartgp( b22bulge, b22d(imax-1),
873 $ rwork(iv2tsn+imax-1-1),
874 $ rwork(iv2tcs+imax-1-1), r )
875 ELSE IF( nu .LT. mu )
THEN
876 CALL dlartgs( b12e(imax-1), b12d(imax), nu,
877 $ rwork(iv2tcs+imax-1-1),
878 $ rwork(iv2tsn+imax-1-1) )
880 CALL dlartgs( b22e(imax-1), b22d(imax), mu,
881 $ rwork(iv2tcs+imax-1-1),
882 $ rwork(iv2tsn+imax-1-1) )
885 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
886 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
887 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
888 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
890 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
891 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
892 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
893 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
900 CALL zlasr(
'R',
'V',
'F', p, imax-imin+1,
901 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
904 CALL zlasr(
'L',
'V',
'F', imax-imin+1, p,
905 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
911 CALL zlasr(
'R',
'V',
'F', m-p, imax-imin+1,
912 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
915 CALL zlasr(
'L',
'V',
'F', imax-imin+1, m-p,
916 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
922 CALL zlasr(
'L',
'V',
'F', imax-imin+1, q,
923 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
924 $ v1t(imin,1), ldv1t )
926 CALL zlasr(
'R',
'V',
'F', q, imax-imin+1,
927 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
928 $ v1t(1,imin), ldv1t )
933 CALL zlasr(
'L',
'V',
'F', imax-imin+1, m-q,
934 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
935 $ v2t(imin,1), ldv2t )
937 CALL zlasr(
'R',
'V',
'F', m-q, imax-imin+1,
938 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
939 $ v2t(1,imin), ldv2t )
945 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 )
THEN
946 b11d(imax) = -b11d(imax)
947 b21d(imax) = -b21d(imax)
950 CALL zscal( q, negonecomplex, v1t(imax,1), ldv1t )
952 CALL zscal( q, negonecomplex, v1t(1,imax), 1 )
959 x1 = cos(phi(imax-1))*b11d(imax) +
960 $ sin(phi(imax-1))*b12e(imax-1)
961 y1 = cos(phi(imax-1))*b21d(imax) +
962 $ sin(phi(imax-1))*b22e(imax-1)
964 theta(imax) = atan2( abs(y1), abs(x1) )
969 IF( b11d(imax)+b12e(imax-1) .LT. 0 )
THEN
970 b12d(imax) = -b12d(imax)
973 CALL zscal( p, negonecomplex, u1(1,imax), 1 )
975 CALL zscal( p, negonecomplex, u1(imax,1), ldu1 )
979 IF( b21d(imax)+b22e(imax-1) .GT. 0 )
THEN
980 b22d(imax) = -b22d(imax)
983 CALL zscal( m-p, negonecomplex, u2(1,imax), 1 )
985 CALL zscal( m-p, negonecomplex, u2(imax,1), ldu2 )
992 IF( b12d(imax)+b22d(imax) .LT. 0 )
THEN
995 CALL zscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
997 CALL zscal( m-q, negonecomplex, v2t(1,imax), 1 )
1005 IF( theta(i) .LT. thresh )
THEN
1007 ELSE IF( theta(i) .GT. piover2-thresh )
THEN
1012 IF( phi(i) .LT. thresh )
THEN
1014 ELSE IF( phi(i) .GT. piover2-thresh )
THEN
1021 IF (imax .GT. 1)
THEN
1022 DO WHILE( phi(imax-1) .EQ. zero )
1024 IF (imax .LE. 1)
EXIT
1027 IF( imin .GT. imax - 1 )
1029 IF (imin .GT. 1)
THEN
1030 DO WHILE (phi(imin-1) .NE. zero)
1032 IF (imin .LE. 1)
EXIT
1047 IF( theta(j) .LT. thetamin )
THEN
1053 IF( mini .NE. i )
THEN
1054 theta(mini) = theta(i)
1058 $
CALL zswap( p, u1(1,i), 1, u1(1,mini), 1 )
1060 $
CALL zswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1062 $
CALL zswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1064 $
CALL zswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1068 $
CALL zswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1070 $
CALL zswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1072 $
CALL zswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1074 $
CALL zswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.
double precision function dlamch(CMACH)
DLAMCH
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine xerbla(SRNAME, INFO)
XERBLA
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 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.
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL