228 INTEGER info, lda, ldv, lwork, m, mv, n, nsweep
229 DOUBLE PRECISION eps, sfmin, tol
233 COMPLEX*16 a( lda, * ), d( n ), v( ldv, * ), work( lwork )
234 DOUBLE PRECISION sva( n )
240 DOUBLE PRECISION zero, half, one
241 parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
242 COMPLEX*16 czero, cone
243 parameter( czero = (0.0d0, 0.0d0), cone = (1.0d0, 0.0d0) )
246 COMPLEX*16 aapq, ompq
247 DOUBLE PRECISION 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, dble, 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(
'ZGSVJ0', -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 =
idamax( n-p+1, sva( p ), 1 ) + p - 1
397 CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
398 IF( rsvec )
CALL zswap( mvl, v( 1, p ), 1,
422 IF( ( sva( p ).LT.rootbig ) .AND.
423 $ ( sva( p ).GT.rootsfmin ) )
THEN
424 sva( p ) =
dznrm2( m, a( 1, p ), 1 )
428 CALL zlassq( 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 = (
zdotc( m, a( 1, p ), 1,
451 $ a( 1, q ), 1 ) / aaqq ) / aapp
453 CALL zcopy( m, a( 1, p ), 1,
455 CALL zlascl(
'G', 0, 0, aapp, one,
456 $ m, 1, work, lda, ierr )
457 aapq =
zdotc( m, work, 1,
458 $ a( 1, q ), 1 ) / aaqq
461 rotok = aapp.LE.( aaqq / small )
462 IF( aapp.GT.( small / aaqq ) )
THEN
463 aapq = (
zdotc( m, a( 1, p ), 1,
464 $ a( 1, q ), 1 ) / aapp ) / aaqq
466 CALL zcopy( m, a( 1, q ), 1,
468 CALL zlascl(
'G', 0, 0, aaqq,
471 aapq =
zdotc( 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 zrot( m, a(1,p), 1, a(1,q), 1,
506 $ cs, conjg(ompq)*t )
508 CALL zrot( 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 zrot( m, a(1,p), 1, a(1,q), 1,
535 $ cs, conjg(ompq)*sn )
537 CALL zrot( mvl, v(1,p), 1,
538 $ v(1,q), 1, cs, conjg(ompq)*sn )
545 CALL zcopy( m, a( 1, p ), 1,
547 CALL zlascl(
'G', 0, 0, aapp, one, m,
550 CALL zlascl(
'G', 0, 0, aaqq, one, m,
551 $ 1, a( 1, q ), lda, ierr )
552 CALL zaxpy( m, -aapq, work, 1,
554 CALL zlascl(
'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 ) =
dznrm2( m, a( 1, q ), 1 )
573 CALL zlassq( 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 =
dznrm2( m, a( 1, p ), 1 )
585 CALL zlassq( 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 = (
zdotc( m, a( 1, p ), 1,
667 $ a( 1, q ), 1 ) / aaqq ) / aapp
669 CALL zcopy( m, a( 1, p ), 1,
671 CALL zlascl(
'G', 0, 0, aapp,
674 aapq =
zdotc( 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 = (
zdotc( m, a( 1, p ), 1,
685 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
688 CALL zcopy( m, a( 1, q ), 1,
690 CALL zlascl(
'G', 0, 0, aaqq,
693 aapq =
zdotc( 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 zrot( m, a(1,p), 1, a(1,q), 1,
722 $ cs, conjg(ompq)*t )
724 CALL zrot( 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 zrot( m, a(1,p), 1, a(1,q), 1,
749 $ cs, conjg(ompq)*sn )
751 CALL zrot( mvl, v(1,p), 1,
752 $ v(1,q), 1, cs, conjg(ompq)*sn )
759 IF( aapp.GT.aaqq )
THEN
760 CALL zcopy( m, a( 1, p ), 1,
762 CALL zlascl(
'G', 0, 0, aapp, one,
765 CALL zlascl(
'G', 0, 0, aaqq, one,
766 $ m, 1, a( 1, q ), lda,
768 CALL zaxpy( m, -aapq, work,
770 CALL zlascl(
'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 zcopy( m, a( 1, q ), 1,
779 CALL zlascl(
'G', 0, 0, aaqq, one,
782 CALL zlascl(
'G', 0, 0, aapp, one,
783 $ m, 1, a( 1, p ), lda,
785 CALL zaxpy( m, -conjg(aapq),
786 $ work, 1, a( 1, p ), 1 )
787 CALL zlascl(
'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 ) =
dznrm2( m, a( 1, q ), 1)
807 CALL zlassq( 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 =
dznrm2( m, a( 1, p ), 1 )
819 CALL zlassq( 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 ) =
dznrm2( m, a( 1, n ), 1 )
885 CALL zlassq( 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( dble( n ) )*
895 $ tol ) .AND. ( dble( n )*mxaapq*mxsinj.LT.tol ) )
THEN
899 IF( notrot.GE.emptsw )
GO TO 1994
918 q =
idamax( n-p+1, sva( p ), 1 ) + p - 1
926 CALL zswap( m, a( 1, p ), 1, a( 1, q ), 1 )
927 IF( rsvec )
CALL zswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
integer function idamax(N, DX, INCX)
IDAMAX
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
complex *16 function zdotc(N, ZX, INCX, ZY, INCY)
ZDOTC
double precision function dznrm2(N, X, INCX)
DZNRM2
subroutine zlassq(N, X, INCX, SCALE, SUMSQ)
ZLASSQ updates a sum of squares represented in scaled form.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY