330 SUBROUTINE cbbcsd( 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 REAL B11D( * ), B11E( * ), B12D( * ), B12E( * ),
346 $ b21d( * ), b21e( * ), b22d( * ), b22e( * ),
347 $ phi( * ), theta( * ), rwork( * )
348 COMPLEX U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
356 parameter ( maxitr = 6 )
357 REAL HUNDRED, MEIGHTH, ONE, PIOVER2, TEN, ZERO
358 parameter ( hundred = 100.0e0, meighth = -0.125e0,
359 $ one = 1.0e0, piover2 = 1.57079632679489662e0,
360 $ ten = 10.0e0, zero = 0.0e0 )
361 COMPLEX NEGONECOMPLEX
362 parameter ( negonecomplex = (-1.0e0,0.0e0) )
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 REAL 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
383 EXTERNAL lsame, slamch
386 INTRINSIC abs, atan2, cos, max, min, sin, sqrt
393 lquery = lrwork .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 lrworkopt = iv2tsn + q - 1
438 lrworkmin = lrworkopt
440 IF( lrwork .LT. lrworkmin .AND. .NOT. lquery )
THEN
445 IF( info .NE. 0 )
THEN
446 CALL xerbla(
'CBBCSD', -info )
448 ELSE IF( lquery )
THEN
454 eps = slamch(
'Epsilon' )
455 unfl = slamch(
'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 slas2( b11d(imax-1), b11e(imax-1), b11d(imax), sigma11,
566 CALL slas2( 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 slartgs( b11d(imin), b11e(imin), mu,
590 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
592 CALL slartgs( b21d(imin), b21e(imin), nu,
593 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1) )
596 temp = rwork(iv1tcs+imin-1)*b11d(imin) +
597 $ rwork(iv1tsn+imin-1)*b11e(imin)
598 b11e(imin) = rwork(iv1tcs+imin-1)*b11e(imin) -
599 $ rwork(iv1tsn+imin-1)*b11d(imin)
601 b11bulge = rwork(iv1tsn+imin-1)*b11d(imin+1)
602 b11d(imin+1) = rwork(iv1tcs+imin-1)*b11d(imin+1)
603 temp = rwork(iv1tcs+imin-1)*b21d(imin) +
604 $ rwork(iv1tsn+imin-1)*b21e(imin)
605 b21e(imin) = rwork(iv1tcs+imin-1)*b21e(imin) -
606 $ rwork(iv1tsn+imin-1)*b21d(imin)
608 b21bulge = rwork(iv1tsn+imin-1)*b21d(imin+1)
609 b21d(imin+1) = rwork(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 slartgp( b11bulge, b11d(imin), rwork(iu1sn+imin-1),
620 $ rwork(iu1cs+imin-1), r )
621 ELSE IF( mu .LE. nu )
THEN
622 CALL slartgs( b11e( imin ), b11d( imin + 1 ), mu,
623 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
625 CALL slartgs( b12d( imin ), b12e( imin ), nu,
626 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1) )
628 IF( b21d(imin)**2+b21bulge**2 .GT. thresh**2 )
THEN
629 CALL slartgp( b21bulge, b21d(imin), rwork(iu2sn+imin-1),
630 $ rwork(iu2cs+imin-1), r )
631 ELSE IF( nu .LT. mu )
THEN
632 CALL slartgs( b21e( imin ), b21d( imin + 1 ), nu,
633 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
635 CALL slartgs( b22d(imin), b22e(imin), mu,
636 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1) )
638 rwork(iu2cs+imin-1) = -rwork(iu2cs+imin-1)
639 rwork(iu2sn+imin-1) = -rwork(iu2sn+imin-1)
641 temp = rwork(iu1cs+imin-1)*b11e(imin) +
642 $ rwork(iu1sn+imin-1)*b11d(imin+1)
643 b11d(imin+1) = rwork(iu1cs+imin-1)*b11d(imin+1) -
644 $ rwork(iu1sn+imin-1)*b11e(imin)
646 IF( imax .GT. imin+1 )
THEN
647 b11bulge = rwork(iu1sn+imin-1)*b11e(imin+1)
648 b11e(imin+1) = rwork(iu1cs+imin-1)*b11e(imin+1)
650 temp = rwork(iu1cs+imin-1)*b12d(imin) +
651 $ rwork(iu1sn+imin-1)*b12e(imin)
652 b12e(imin) = rwork(iu1cs+imin-1)*b12e(imin) -
653 $ rwork(iu1sn+imin-1)*b12d(imin)
655 b12bulge = rwork(iu1sn+imin-1)*b12d(imin+1)
656 b12d(imin+1) = rwork(iu1cs+imin-1)*b12d(imin+1)
657 temp = rwork(iu2cs+imin-1)*b21e(imin) +
658 $ rwork(iu2sn+imin-1)*b21d(imin+1)
659 b21d(imin+1) = rwork(iu2cs+imin-1)*b21d(imin+1) -
660 $ rwork(iu2sn+imin-1)*b21e(imin)
662 IF( imax .GT. imin+1 )
THEN
663 b21bulge = rwork(iu2sn+imin-1)*b21e(imin+1)
664 b21e(imin+1) = rwork(iu2cs+imin-1)*b21e(imin+1)
666 temp = rwork(iu2cs+imin-1)*b22d(imin) +
667 $ rwork(iu2sn+imin-1)*b22e(imin)
668 b22e(imin) = rwork(iu2cs+imin-1)*b22e(imin) -
669 $ rwork(iu2sn+imin-1)*b22d(imin)
671 b22bulge = rwork(iu2sn+imin-1)*b22d(imin+1)
672 b22d(imin+1) = rwork(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 slartgp( x2, x1, rwork(iv1tsn+i-1),
703 $ rwork(iv1tcs+i-1), r )
704 ELSE IF( .NOT. restart11 .AND. restart21 )
THEN
705 CALL slartgp( b11bulge, b11e(i-1), rwork(iv1tsn+i-1),
706 $ rwork(iv1tcs+i-1), r )
707 ELSE IF( restart11 .AND. .NOT. restart21 )
THEN
708 CALL slartgp( b21bulge, b21e(i-1), rwork(iv1tsn+i-1),
709 $ rwork(iv1tcs+i-1), r )
710 ELSE IF( mu .LE. nu )
THEN
711 CALL slartgs( b11d(i), b11e(i), mu, rwork(iv1tcs+i-1),
712 $ rwork(iv1tsn+i-1) )
714 CALL slartgs( b21d(i), b21e(i), nu, rwork(iv1tcs+i-1),
715 $ rwork(iv1tsn+i-1) )
717 rwork(iv1tcs+i-1) = -rwork(iv1tcs+i-1)
718 rwork(iv1tsn+i-1) = -rwork(iv1tsn+i-1)
719 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
720 CALL slartgp( y2, y1, rwork(iv2tsn+i-1-1),
721 $ rwork(iv2tcs+i-1-1), r )
722 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
723 CALL slartgp( b12bulge, b12d(i-1), rwork(iv2tsn+i-1-1),
724 $ rwork(iv2tcs+i-1-1), r )
725 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
726 CALL slartgp( b22bulge, b22d(i-1), rwork(iv2tsn+i-1-1),
727 $ rwork(iv2tcs+i-1-1), r )
728 ELSE IF( nu .LT. mu )
THEN
729 CALL slartgs( b12e(i-1), b12d(i), nu,
730 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
732 CALL slartgs( b22e(i-1), b22d(i), mu,
733 $ rwork(iv2tcs+i-1-1), rwork(iv2tsn+i-1-1) )
736 temp = rwork(iv1tcs+i-1)*b11d(i) + rwork(iv1tsn+i-1)*b11e(i)
737 b11e(i) = rwork(iv1tcs+i-1)*b11e(i) -
738 $ rwork(iv1tsn+i-1)*b11d(i)
740 b11bulge = rwork(iv1tsn+i-1)*b11d(i+1)
741 b11d(i+1) = rwork(iv1tcs+i-1)*b11d(i+1)
742 temp = rwork(iv1tcs+i-1)*b21d(i) + rwork(iv1tsn+i-1)*b21e(i)
743 b21e(i) = rwork(iv1tcs+i-1)*b21e(i) -
744 $ rwork(iv1tsn+i-1)*b21d(i)
746 b21bulge = rwork(iv1tsn+i-1)*b21d(i+1)
747 b21d(i+1) = rwork(iv1tcs+i-1)*b21d(i+1)
748 temp = rwork(iv2tcs+i-1-1)*b12e(i-1) +
749 $ rwork(iv2tsn+i-1-1)*b12d(i)
750 b12d(i) = rwork(iv2tcs+i-1-1)*b12d(i) -
751 $ rwork(iv2tsn+i-1-1)*b12e(i-1)
753 b12bulge = rwork(iv2tsn+i-1-1)*b12e(i)
754 b12e(i) = rwork(iv2tcs+i-1-1)*b12e(i)
755 temp = rwork(iv2tcs+i-1-1)*b22e(i-1) +
756 $ rwork(iv2tsn+i-1-1)*b22d(i)
757 b22d(i) = rwork(iv2tcs+i-1-1)*b22d(i) -
758 $ rwork(iv2tsn+i-1-1)*b22e(i-1)
760 b22bulge = rwork(iv2tsn+i-1-1)*b22e(i)
761 b22e(i) = rwork(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 slartgp( x2, x1, rwork(iu1sn+i-1), rwork(iu1cs+i-1),
787 ELSE IF( .NOT. restart11 .AND. restart12 )
THEN
788 CALL slartgp( b11bulge, b11d(i), rwork(iu1sn+i-1),
789 $ rwork(iu1cs+i-1), r )
790 ELSE IF( restart11 .AND. .NOT. restart12 )
THEN
791 CALL slartgp( b12bulge, b12e(i-1), rwork(iu1sn+i-1),
792 $ rwork(iu1cs+i-1), r )
793 ELSE IF( mu .LE. nu )
THEN
794 CALL slartgs( b11e(i), b11d(i+1), mu, rwork(iu1cs+i-1),
797 CALL slartgs( b12d(i), b12e(i), nu, rwork(iu1cs+i-1),
800 IF( .NOT. restart21 .AND. .NOT. restart22 )
THEN
801 CALL slartgp( y2, y1, rwork(iu2sn+i-1), rwork(iu2cs+i-1),
803 ELSE IF( .NOT. restart21 .AND. restart22 )
THEN
804 CALL slartgp( b21bulge, b21d(i), rwork(iu2sn+i-1),
805 $ rwork(iu2cs+i-1), r )
806 ELSE IF( restart21 .AND. .NOT. restart22 )
THEN
807 CALL slartgp( b22bulge, b22e(i-1), rwork(iu2sn+i-1),
808 $ rwork(iu2cs+i-1), r )
809 ELSE IF( nu .LT. mu )
THEN
810 CALL slartgs( b21e(i), b21e(i+1), nu, rwork(iu2cs+i-1),
813 CALL slartgs( b22d(i), b22e(i), mu, rwork(iu2cs+i-1),
816 rwork(iu2cs+i-1) = -rwork(iu2cs+i-1)
817 rwork(iu2sn+i-1) = -rwork(iu2sn+i-1)
819 temp = rwork(iu1cs+i-1)*b11e(i) + rwork(iu1sn+i-1)*b11d(i+1)
820 b11d(i+1) = rwork(iu1cs+i-1)*b11d(i+1) -
821 $ rwork(iu1sn+i-1)*b11e(i)
823 IF( i .LT. imax - 1 )
THEN
824 b11bulge = rwork(iu1sn+i-1)*b11e(i+1)
825 b11e(i+1) = rwork(iu1cs+i-1)*b11e(i+1)
827 temp = rwork(iu2cs+i-1)*b21e(i) + rwork(iu2sn+i-1)*b21d(i+1)
828 b21d(i+1) = rwork(iu2cs+i-1)*b21d(i+1) -
829 $ rwork(iu2sn+i-1)*b21e(i)
831 IF( i .LT. imax - 1 )
THEN
832 b21bulge = rwork(iu2sn+i-1)*b21e(i+1)
833 b21e(i+1) = rwork(iu2cs+i-1)*b21e(i+1)
835 temp = rwork(iu1cs+i-1)*b12d(i) + rwork(iu1sn+i-1)*b12e(i)
836 b12e(i) = rwork(iu1cs+i-1)*b12e(i) -
837 $ rwork(iu1sn+i-1)*b12d(i)
839 b12bulge = rwork(iu1sn+i-1)*b12d(i+1)
840 b12d(i+1) = rwork(iu1cs+i-1)*b12d(i+1)
841 temp = rwork(iu2cs+i-1)*b22d(i) + rwork(iu2sn+i-1)*b22e(i)
842 b22e(i) = rwork(iu2cs+i-1)*b22e(i) -
843 $ rwork(iu2sn+i-1)*b22d(i)
845 b22bulge = rwork(iu2sn+i-1)*b22d(i+1)
846 b22d(i+1) = rwork(iu2cs+i-1)*b22d(i+1)
852 x1 = sin(theta(imax-1))*b11e(imax-1) +
853 $ cos(theta(imax-1))*b21e(imax-1)
854 y1 = sin(theta(imax-1))*b12d(imax-1) +
855 $ cos(theta(imax-1))*b22d(imax-1)
856 y2 = sin(theta(imax-1))*b12bulge + cos(theta(imax-1))*b22bulge
858 phi(imax-1) = atan2( abs(x1), sqrt(y1**2+y2**2) )
862 restart12 = b12d(imax-1)**2 + b12bulge**2 .LE. thresh**2
863 restart22 = b22d(imax-1)**2 + b22bulge**2 .LE. thresh**2
865 IF( .NOT. restart12 .AND. .NOT. restart22 )
THEN
866 CALL slartgp( y2, y1, rwork(iv2tsn+imax-1-1),
867 $ rwork(iv2tcs+imax-1-1), r )
868 ELSE IF( .NOT. restart12 .AND. restart22 )
THEN
869 CALL slartgp( b12bulge, b12d(imax-1),
870 $ rwork(iv2tsn+imax-1-1),
871 $ rwork(iv2tcs+imax-1-1), r )
872 ELSE IF( restart12 .AND. .NOT. restart22 )
THEN
873 CALL slartgp( b22bulge, b22d(imax-1),
874 $ rwork(iv2tsn+imax-1-1),
875 $ rwork(iv2tcs+imax-1-1), r )
876 ELSE IF( nu .LT. mu )
THEN
877 CALL slartgs( b12e(imax-1), b12d(imax), nu,
878 $ rwork(iv2tcs+imax-1-1),
879 $ rwork(iv2tsn+imax-1-1) )
881 CALL slartgs( b22e(imax-1), b22d(imax), mu,
882 $ rwork(iv2tcs+imax-1-1),
883 $ rwork(iv2tsn+imax-1-1) )
886 temp = rwork(iv2tcs+imax-1-1)*b12e(imax-1) +
887 $ rwork(iv2tsn+imax-1-1)*b12d(imax)
888 b12d(imax) = rwork(iv2tcs+imax-1-1)*b12d(imax) -
889 $ rwork(iv2tsn+imax-1-1)*b12e(imax-1)
891 temp = rwork(iv2tcs+imax-1-1)*b22e(imax-1) +
892 $ rwork(iv2tsn+imax-1-1)*b22d(imax)
893 b22d(imax) = rwork(iv2tcs+imax-1-1)*b22d(imax) -
894 $ rwork(iv2tsn+imax-1-1)*b22e(imax-1)
901 CALL clasr(
'R',
'V',
'F', p, imax-imin+1,
902 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
905 CALL clasr(
'L',
'V',
'F', imax-imin+1, p,
906 $ rwork(iu1cs+imin-1), rwork(iu1sn+imin-1),
912 CALL clasr(
'R',
'V',
'F', m-p, imax-imin+1,
913 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
916 CALL clasr(
'L',
'V',
'F', imax-imin+1, m-p,
917 $ rwork(iu2cs+imin-1), rwork(iu2sn+imin-1),
923 CALL clasr(
'L',
'V',
'F', imax-imin+1, q,
924 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
925 $ v1t(imin,1), ldv1t )
927 CALL clasr(
'R',
'V',
'F', q, imax-imin+1,
928 $ rwork(iv1tcs+imin-1), rwork(iv1tsn+imin-1),
929 $ v1t(1,imin), ldv1t )
934 CALL clasr(
'L',
'V',
'F', imax-imin+1, m-q,
935 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
936 $ v2t(imin,1), ldv2t )
938 CALL clasr(
'R',
'V',
'F', m-q, imax-imin+1,
939 $ rwork(iv2tcs+imin-1), rwork(iv2tsn+imin-1),
940 $ v2t(1,imin), ldv2t )
946 IF( b11e(imax-1)+b21e(imax-1) .GT. 0 )
THEN
947 b11d(imax) = -b11d(imax)
948 b21d(imax) = -b21d(imax)
951 CALL cscal( q, negonecomplex, v1t(imax,1), ldv1t )
953 CALL cscal( q, negonecomplex, v1t(1,imax), 1 )
960 x1 = cos(phi(imax-1))*b11d(imax) +
961 $ sin(phi(imax-1))*b12e(imax-1)
962 y1 = cos(phi(imax-1))*b21d(imax) +
963 $ sin(phi(imax-1))*b22e(imax-1)
965 theta(imax) = atan2( abs(y1), abs(x1) )
970 IF( b11d(imax)+b12e(imax-1) .LT. 0 )
THEN
971 b12d(imax) = -b12d(imax)
974 CALL cscal( p, negonecomplex, u1(1,imax), 1 )
976 CALL cscal( p, negonecomplex, u1(imax,1), ldu1 )
980 IF( b21d(imax)+b22e(imax-1) .GT. 0 )
THEN
981 b22d(imax) = -b22d(imax)
984 CALL cscal( m-p, negonecomplex, u2(1,imax), 1 )
986 CALL cscal( m-p, negonecomplex, u2(imax,1), ldu2 )
993 IF( b12d(imax)+b22d(imax) .LT. 0 )
THEN
996 CALL cscal( m-q, negonecomplex, v2t(imax,1), ldv2t )
998 CALL cscal( m-q, negonecomplex, v2t(1,imax), 1 )
1006 IF( theta(i) .LT. thresh )
THEN
1008 ELSE IF( theta(i) .GT. piover2-thresh )
THEN
1013 IF( phi(i) .LT. thresh )
THEN
1015 ELSE IF( phi(i) .GT. piover2-thresh )
THEN
1022 IF (imax .GT. 1)
THEN
1023 DO WHILE( phi(imax-1) .EQ. zero )
1025 IF (imax .LE. 1)
EXIT
1028 IF( imin .GT. imax - 1 )
1030 IF (imin .GT. 1)
THEN
1031 DO WHILE (phi(imin-1) .NE. zero)
1033 IF (imin .LE. 1)
EXIT
1048 IF( theta(j) .LT. thetamin )
THEN
1054 IF( mini .NE. i )
THEN
1055 theta(mini) = theta(i)
1059 $
CALL cswap( p, u1(1,i), 1, u1(1,mini), 1 )
1061 $
CALL cswap( m-p, u2(1,i), 1, u2(1,mini), 1 )
1063 $
CALL cswap( q, v1t(i,1), ldv1t, v1t(mini,1), ldv1t )
1065 $
CALL cswap( m-q, v2t(i,1), ldv2t, v2t(mini,1),
1069 $
CALL cswap( p, u1(i,1), ldu1, u1(mini,1), ldu1 )
1071 $
CALL cswap( m-p, u2(i,1), ldu2, u2(mini,1), ldu2 )
1073 $
CALL cswap( q, v1t(1,i), 1, v1t(1,mini), 1 )
1075 $
CALL cswap( m-q, v2t(1,i), 1, v2t(1,mini), 1 )
subroutine cbbcsd(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)
CBBCSD
subroutine slartgs(X, Y, SIGMA, CS, SN)
SLARTGS generates a plane rotation designed to introduce a bulge in implicit QR iteration for the bid...
subroutine slas2(F, G, H, SSMIN, SSMAX)
SLAS2 computes singular values of a 2-by-2 triangular matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cscal(N, CA, CX, INCX)
CSCAL
subroutine clasr(SIDE, PIVOT, DIRECT, M, N, C, S, A, LDA)
CLASR applies a sequence of plane rotations to a general rectangular matrix.
subroutine slartgp(F, G, CS, SN, R)
SLARTGP generates a plane rotation so that the diagonal is nonnegative.
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP