228 INTEGER info, lda, ldv, lwork, m, mv, n, nsweep
233 COMPLEX a( lda, * ), d( n ), v( ldv, * ), work( lwork )
241 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
243 parameter( czero = (0.0e0, 0.0e0), cone = (1.0e0, 0.0e0) )
247 REAL aapp, aapp0, aapq1, aaqq, apoaq, aqoap, big,
248 $ bigtheta, cs, mxaapq, mxsinj, rootbig, rooteps,
249 $ rootsfmin, roottol, small, sn, t, temp1, theta,
251 INTEGER blskip, emptsw, i, ibr, ierr, igl, ijblsk, ir1,
252 $ iswrot, jbc, jgl, kbl, lkahead, mvl, nbl,
253 $ notrot, p, pskipped, q, rowskip, swband
254 LOGICAL applv, rotok, rsvec
258 INTRINSIC abs, max, conjg,
REAL, min, sign, sqrt
279 applv =
lsame( jobv,
'A' )
280 rsvec =
lsame( jobv,
'V' )
281 IF( .NOT.( rsvec .OR. applv .OR.
lsame( jobv,
'N' ) ) )
THEN
283 ELSE IF( m.LT.0 )
THEN
285 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
287 ELSE IF( lda.LT.m )
THEN
289 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
291 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
292 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
294 ELSE IF( tol.LE.eps )
THEN
296 ELSE IF( nsweep.LT.0 )
THEN
298 ELSE IF( lwork.LT.m )
THEN
306 CALL xerbla(
'CGSVJ0', -info )
312 ELSE IF( applv )
THEN
315 rsvec = rsvec .OR. applv
317 rooteps = sqrt( eps )
318 rootsfmin = sqrt( sfmin )
321 rootbig = one / rootsfmin
322 bigtheta = one / rooteps
323 roottol = sqrt( tol )
327 emptsw = ( n*( n-1 ) ) / 2
348 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
353 rowskip = min( 5, kbl )
367 DO 1993 i = 1, nsweep
385 igl = ( ibr-1 )*kbl + 1
387 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
391 DO 2001 p = igl, min( igl+kbl-1, n-1 )
395 q =
isamax( n-p+1, sva( p ), 1 ) + p - 1
397 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
398 IF( rsvec )
CALL cswap( mvl, v( 1, p ), 1,
422 IF( ( sva( p ).LT.rootbig ) .AND.
423 $ ( sva( p ).GT.rootsfmin ) )
THEN
424 sva( p ) =
scnrm2( m, a( 1, p ), 1 )
428 CALL classq( m, a( 1, p ), 1, temp1, aapp )
429 sva( p ) = temp1*sqrt( aapp )
436 IF( aapp.GT.zero )
THEN
440 DO 2002 q = p + 1, min( igl+kbl-1, n )
444 IF( aaqq.GT.zero )
THEN
447 IF( aaqq.GE.one )
THEN
448 rotok = ( small*aapp ).LE.aaqq
449 IF( aapp.LT.( big / aaqq ) )
THEN
450 aapq = (
cdotc( m, a( 1, p ), 1,
451 $ a( 1, q ), 1 ) / aaqq ) / aapp
453 CALL ccopy( m, a( 1, p ), 1,
455 CALL clascl(
'G', 0, 0, aapp, one,
456 $ m, 1, work, lda, ierr )
457 aapq =
cdotc( m, work, 1,
458 $ a( 1, q ), 1 ) / aaqq
461 rotok = aapp.LE.( aaqq / small )
462 IF( aapp.GT.( small / aaqq ) )
THEN
463 aapq = (
cdotc( m, a( 1, p ), 1,
464 $ a( 1, q ), 1 ) / aapp ) / aaqq
466 CALL ccopy( m, a( 1, q ), 1,
468 CALL clascl(
'G', 0, 0, aaqq,
471 aapq =
cdotc( m, a( 1, p ), 1,
478 mxaapq = max( mxaapq, -aapq1 )
482 IF( abs( aapq1 ).GT.tol )
THEN
483 ompq = aapq / abs(aapq)
498 theta = -half*abs( aqoap-apoaq )/aapq1
500 IF( abs( theta ).GT.bigtheta )
THEN
505 CALL crot( m, a(1,p), 1, a(1,q), 1,
506 $ cs, conjg(ompq)*t )
508 CALL crot( mvl, v(1,p), 1,
509 $ v(1,q), 1, cs, conjg(ompq)*t )
512 sva( q ) = aaqq*sqrt( max( zero,
513 $ one+t*apoaq*aapq1 ) )
514 aapp = aapp*sqrt( max( zero,
515 $ one-t*aqoap*aapq1 ) )
516 mxsinj = max( mxsinj, abs( t ) )
522 thsign = -sign( one, aapq1 )
523 t = one / ( theta+thsign*
524 $ sqrt( one+theta*theta ) )
525 cs = sqrt( one / ( one+t*t ) )
528 mxsinj = max( mxsinj, abs( sn ) )
529 sva( q ) = aaqq*sqrt( max( zero,
530 $ one+t*apoaq*aapq1 ) )
531 aapp = aapp*sqrt( max( zero,
532 $ one-t*aqoap*aapq1 ) )
534 CALL crot( m, a(1,p), 1, a(1,q), 1,
535 $ cs, conjg(ompq)*sn )
537 CALL crot( mvl, v(1,p), 1,
538 $ v(1,q), 1, cs, conjg(ompq)*sn )
545 CALL ccopy( m, a( 1, p ), 1,
547 CALL clascl(
'G', 0, 0, aapp, one, m,
550 CALL clascl(
'G', 0, 0, aaqq, one, m,
551 $ 1, a( 1, q ), lda, ierr )
552 CALL caxpy( m, -aapq, work, 1,
554 CALL clascl(
'G', 0, 0, one, aaqq, m,
555 $ 1, a( 1, q ), lda, ierr )
556 sva( q ) = aaqq*sqrt( max( zero,
557 $ one-aapq1*aapq1 ) )
558 mxsinj = max( mxsinj, sfmin )
565 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
567 IF( ( aaqq.LT.rootbig ) .AND.
568 $ ( aaqq.GT.rootsfmin ) )
THEN
569 sva( q ) =
scnrm2( m, a( 1, q ), 1 )
573 CALL classq( m, a( 1, q ), 1, t,
575 sva( q ) = t*sqrt( aaqq )
578 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
579 IF( ( aapp.LT.rootbig ) .AND.
580 $ ( aapp.GT.rootsfmin ) )
THEN
581 aapp =
scnrm2( m, a( 1, p ), 1 )
585 CALL classq( m, a( 1, p ), 1, t,
587 aapp = t*sqrt( aapp )
594 IF( ir1.EQ.0 )notrot = notrot + 1
596 pskipped = pskipped + 1
600 IF( ir1.EQ.0 )notrot = notrot + 1
601 pskipped = pskipped + 1
604 IF( ( i.LE.swband ) .AND.
605 $ ( pskipped.GT.rowskip ) )
THEN
606 IF( ir1.EQ.0 )aapp = -aapp
621 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
622 $ notrot = notrot + min( igl+kbl-1, n ) - p
633 igl = ( ibr-1 )*kbl + 1
635 DO 2010 jbc = ibr + 1, nbl
637 jgl = ( jbc-1 )*kbl + 1
642 DO 2100 p = igl, min( igl+kbl-1, n )
645 IF( aapp.GT.zero )
THEN
649 DO 2200 q = jgl, min( jgl+kbl-1, n )
652 IF( aaqq.GT.zero )
THEN
659 IF( aaqq.GE.one )
THEN
660 IF( aapp.GE.aaqq )
THEN
661 rotok = ( small*aapp ).LE.aaqq
663 rotok = ( small*aaqq ).LE.aapp
665 IF( aapp.LT.( big / aaqq ) )
THEN
666 aapq = (
cdotc( m, a( 1, p ), 1,
667 $ a( 1, q ), 1 ) / aaqq ) / aapp
669 CALL ccopy( m, a( 1, p ), 1,
671 CALL clascl(
'G', 0, 0, aapp,
674 aapq =
cdotc( m, work, 1,
675 $ a( 1, q ), 1 ) / aaqq
678 IF( aapp.GE.aaqq )
THEN
679 rotok = aapp.LE.( aaqq / small )
681 rotok = aaqq.LE.( aapp / small )
683 IF( aapp.GT.( small / aaqq ) )
THEN
684 aapq = (
cdotc( m, a( 1, p ), 1,
685 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
688 CALL ccopy( m, a( 1, q ), 1,
690 CALL clascl(
'G', 0, 0, aaqq,
693 aapq =
cdotc( m, a( 1, p ), 1,
700 mxaapq = max( mxaapq, -aapq1 )
704 IF( abs( aapq1 ).GT.tol )
THEN
705 ompq = aapq / abs(aapq)
715 theta = -half*abs( aqoap-apoaq )/ aapq1
716 IF( aaqq.GT.aapp0 )theta = -theta
718 IF( abs( theta ).GT.bigtheta )
THEN
721 CALL crot( m, a(1,p), 1, a(1,q), 1,
722 $ cs, conjg(ompq)*t )
724 CALL crot( mvl, v(1,p), 1,
725 $ v(1,q), 1, cs, conjg(ompq)*t )
727 sva( q ) = aaqq*sqrt( max( zero,
728 $ one+t*apoaq*aapq1 ) )
729 aapp = aapp*sqrt( max( zero,
730 $ one-t*aqoap*aapq1 ) )
731 mxsinj = max( mxsinj, abs( t ) )
736 thsign = -sign( one, aapq1 )
737 IF( aaqq.GT.aapp0 )thsign = -thsign
738 t = one / ( theta+thsign*
739 $ sqrt( one+theta*theta ) )
740 cs = sqrt( one / ( one+t*t ) )
742 mxsinj = max( mxsinj, abs( sn ) )
743 sva( q ) = aaqq*sqrt( max( zero,
744 $ one+t*apoaq*aapq1 ) )
745 aapp = aapp*sqrt( max( zero,
746 $ one-t*aqoap*aapq1 ) )
748 CALL crot( m, a(1,p), 1, a(1,q), 1,
749 $ cs, conjg(ompq)*sn )
751 CALL crot( mvl, v(1,p), 1,
752 $ v(1,q), 1, cs, conjg(ompq)*sn )
759 IF( aapp.GT.aaqq )
THEN
760 CALL ccopy( m, a( 1, p ), 1,
762 CALL clascl(
'G', 0, 0, aapp, one,
765 CALL clascl(
'G', 0, 0, aaqq, one,
766 $ m, 1, a( 1, q ), lda,
768 CALL caxpy( m, -aapq, work,
770 CALL clascl(
'G', 0, 0, one, aaqq,
771 $ m, 1, a( 1, q ), lda,
773 sva( q ) = aaqq*sqrt( max( zero,
774 $ one-aapq1*aapq1 ) )
775 mxsinj = max( mxsinj, sfmin )
777 CALL ccopy( m, a( 1, q ), 1,
779 CALL clascl(
'G', 0, 0, aaqq, one,
782 CALL clascl(
'G', 0, 0, aapp, one,
783 $ m, 1, a( 1, p ), lda,
785 CALL caxpy( m, -conjg(aapq),
786 $ work, 1, a( 1, p ), 1 )
787 CALL clascl(
'G', 0, 0, one, aapp,
788 $ m, 1, a( 1, p ), lda,
790 sva( p ) = aapp*sqrt( max( zero,
791 $ one-aapq1*aapq1 ) )
792 mxsinj = max( mxsinj, sfmin )
799 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
801 IF( ( aaqq.LT.rootbig ) .AND.
802 $ ( aaqq.GT.rootsfmin ) )
THEN
803 sva( q ) =
scnrm2( m, a( 1, q ), 1)
807 CALL classq( m, a( 1, q ), 1, t,
809 sva( q ) = t*sqrt( aaqq )
812 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
813 IF( ( aapp.LT.rootbig ) .AND.
814 $ ( aapp.GT.rootsfmin ) )
THEN
815 aapp =
scnrm2( m, a( 1, p ), 1 )
819 CALL classq( m, a( 1, p ), 1, t,
821 aapp = t*sqrt( aapp )
829 pskipped = pskipped + 1
834 pskipped = pskipped + 1
838 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
844 IF( ( i.LE.swband ) .AND.
845 $ ( pskipped.GT.rowskip ) )
THEN
859 IF( aapp.EQ.zero )notrot = notrot +
860 $ min( jgl+kbl-1, n ) - jgl + 1
861 IF( aapp.LT.zero )notrot = 0
871 DO 2012 p = igl, min( igl+kbl-1, n )
872 sva( p ) = abs( sva( p ) )
879 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
881 sva( n ) =
scnrm2( m, a( 1, n ), 1 )
885 CALL classq( m, a( 1, n ), 1, t, aapp )
886 sva( n ) = t*sqrt( aapp )
891 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
892 $ ( iswrot.LE.n ) ) )swband = i
894 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt(
REAL( N ) )*
895 $ tol ) .AND. (
REAL( n )*mxaapq*mxsinj.LT.tol ) ) then
899 IF( notrot.GE.emptsw )
GO TO 1994
918 q =
isamax( n-p+1, sva( p ), 1 ) + p - 1
926 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
927 IF( rsvec )
CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
real function scnrm2(N, X, INCX)
SCNRM2
subroutine classq(N, X, INCX, SCALE, SUMSQ)
CLASSQ updates a sum of squares represented in scaled form.
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
integer function isamax(N, SX, INCX)
ISAMAX
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY