330 SUBROUTINE zbbcsd( 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, rwork, lrwork, info )
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
380 DOUBLE PRECISION DLAMCH
382 EXTERNAL lsame, dlamch
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
453 eps = dlamch(
'Epsilon' )
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 zbbcsd(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, RWORK, LRWORK, INFO)
ZBBCSD
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 zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlartgp(F, G, CS, SN, R)
DLARTGP generates a plane rotation so that the diagonal is nonnegative.
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
subroutine zlasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
ZLASR applies a sequence of plane rotations to a general rectangular matrix.