227 INTEGER info, lda, ldv, lwork, m, mv, n, nsweep
228 DOUBLE PRECISION eps, sfmin, tol
232 DOUBLE PRECISION a( lda, * ), sva( n ), d( n ), v( ldv, * ),
239 DOUBLE PRECISION zero, half, one
240 parameter ( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
243 DOUBLE PRECISION aapp, aapp0, aapq, aaqq, apoaq, aqoap, big,
244 $ bigtheta, cs, mxaapq, mxsinj, rootbig, rooteps,
245 $ rootsfmin, roottol, small, sn, t, temp1, theta,
247 INTEGER blskip, emptsw, i, ibr, ierr, igl, ijblsk, ir1,
248 $ iswrot, jbc, jgl, kbl, lkahead, mvl, nbl,
249 $ notrot, p, pskipped, q, rowskip, swband
250 LOGICAL applv, rotok, rsvec
253 DOUBLE PRECISION fastr( 5 )
256 INTRINSIC dabs, max, dble, min, dsign, dsqrt
271 applv =
lsame( jobv,
'A' )
272 rsvec =
lsame( jobv,
'V' )
273 IF( .NOT.( rsvec .OR. applv .OR.
lsame( jobv,
'N' ) ) )
THEN
275 ELSE IF( m.LT.0 )
THEN
277 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
279 ELSE IF( lda.LT.m )
THEN
281 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
283 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
284 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
286 ELSE IF( tol.LE.eps )
THEN
288 ELSE IF( nsweep.LT.0 )
THEN
290 ELSE IF( lwork.LT.m )
THEN
298 CALL xerbla(
'DGSVJ0', -info )
304 ELSE IF( applv )
THEN
307 rsvec = rsvec .OR. applv
309 rooteps = dsqrt( eps )
310 rootsfmin = dsqrt( sfmin )
313 rootbig = one / rootsfmin
314 bigtheta = one / rooteps
315 roottol = dsqrt( tol )
319 emptsw = ( n*( n-1 ) ) / 2
339 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
341 blskip = ( kbl**2 ) + 1
344 rowskip = min( 5, kbl )
352 DO 1993 i = 1, nsweep
364 igl = ( ibr-1 )*kbl + 1
366 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
370 DO 2001 p = igl, min( igl+kbl-1, n-1 )
373 q =
idamax( n-p+1, sva( p ), 1 ) + p - 1
375 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
376 IF( rsvec )
CALL dswap( mvl, v( 1, p ), 1,
400 IF( ( sva( p ).LT.rootbig ) .AND.
401 $ ( sva( p ).GT.rootsfmin ) )
THEN
402 sva( p ) =
dnrm2( m, a( 1, p ), 1 )*d( p )
406 CALL dlassq( m, a( 1, p ), 1, temp1, aapp )
407 sva( p ) = temp1*dsqrt( aapp )*d( p )
415 IF( aapp.GT.zero )
THEN
419 DO 2002 q = p + 1, min( igl+kbl-1, n )
423 IF( aaqq.GT.zero )
THEN
426 IF( aaqq.GE.one )
THEN
427 rotok = ( small*aapp ).LE.aaqq
428 IF( aapp.LT.( big / aaqq ) )
THEN
429 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
430 $ q ), 1 )*d( p )*d( q ) / aaqq )
433 CALL dcopy( m, a( 1, p ), 1, work, 1 )
434 CALL dlascl(
'G', 0, 0, aapp, d( p ),
435 $ m, 1, work, lda, ierr )
436 aapq =
ddot( m, work, 1, a( 1, q ),
440 rotok = aapp.LE.( aaqq / small )
441 IF( aapp.GT.( small / aaqq ) )
THEN
442 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
443 $ q ), 1 )*d( p )*d( q ) / aaqq )
446 CALL dcopy( m, a( 1, q ), 1, work, 1 )
447 CALL dlascl(
'G', 0, 0, aaqq, d( q ),
448 $ m, 1, work, lda, ierr )
449 aapq =
ddot( m, work, 1, a( 1, p ),
454 mxaapq = max( mxaapq, dabs( aapq ) )
458 IF( dabs( aapq ).GT.tol )
THEN
473 theta = -half*dabs( aqoap-apoaq )/aapq
475 IF( dabs( theta ).GT.bigtheta )
THEN
478 fastr( 3 ) = t*d( p ) / d( q )
479 fastr( 4 ) = -t*d( q ) / d( p )
480 CALL drotm( m, a( 1, p ), 1,
481 $ a( 1, q ), 1, fastr )
482 IF( rsvec )
CALL drotm( mvl,
486 sva( q ) = aaqq*dsqrt( max( zero,
487 $ one+t*apoaq*aapq ) )
488 aapp = aapp*dsqrt( max( zero,
489 $ one-t*aqoap*aapq ) )
490 mxsinj = max( mxsinj, dabs( t ) )
496 thsign = -dsign( one, aapq )
497 t = one / ( theta+thsign*
498 $ dsqrt( one+theta*theta ) )
499 cs = dsqrt( one / ( one+t*t ) )
502 mxsinj = max( mxsinj, dabs( sn ) )
503 sva( q ) = aaqq*dsqrt( max( zero,
504 $ one+t*apoaq*aapq ) )
505 aapp = aapp*dsqrt( max( zero,
506 $ one-t*aqoap*aapq ) )
508 apoaq = d( p ) / d( q )
509 aqoap = d( q ) / d( p )
510 IF( d( p ).GE.one )
THEN
511 IF( d( q ).GE.one )
THEN
513 fastr( 4 ) = -t*aqoap
516 CALL drotm( m, a( 1, p ), 1,
519 IF( rsvec )
CALL drotm( mvl,
520 $ v( 1, p ), 1, v( 1, q ),
523 CALL daxpy( m, -t*aqoap,
526 CALL daxpy( m, cs*sn*apoaq,
532 CALL daxpy( mvl, -t*aqoap,
542 IF( d( q ).GE.one )
THEN
543 CALL daxpy( m, t*apoaq,
546 CALL daxpy( m, -cs*sn*aqoap,
552 CALL daxpy( mvl, t*apoaq,
561 IF( d( p ).GE.d( q ) )
THEN
562 CALL daxpy( m, -t*aqoap,
565 CALL daxpy( m, cs*sn*apoaq,
581 CALL daxpy( m, t*apoaq,
592 $ t*apoaq, v( 1, p ),
606 CALL dcopy( m, a( 1, p ), 1, work, 1 )
607 CALL dlascl(
'G', 0, 0, aapp, one, m,
608 $ 1, work, lda, ierr )
609 CALL dlascl(
'G', 0, 0, aaqq, one, m,
610 $ 1, a( 1, q ), lda, ierr )
611 temp1 = -aapq*d( p ) / d( q )
612 CALL daxpy( m, temp1, work, 1,
614 CALL dlascl(
'G', 0, 0, one, aaqq, m,
615 $ 1, a( 1, q ), lda, ierr )
616 sva( q ) = aaqq*dsqrt( max( zero,
618 mxsinj = max( mxsinj, sfmin )
624 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
626 IF( ( aaqq.LT.rootbig ) .AND.
627 $ ( aaqq.GT.rootsfmin ) )
THEN
628 sva( q ) =
dnrm2( m, a( 1, q ), 1 )*
633 CALL dlassq( m, a( 1, q ), 1, t,
635 sva( q ) = t*dsqrt( aaqq )*d( q )
638 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
639 IF( ( aapp.LT.rootbig ) .AND.
640 $ ( aapp.GT.rootsfmin ) )
THEN
641 aapp =
dnrm2( m, a( 1, p ), 1 )*
646 CALL dlassq( m, a( 1, p ), 1, t,
648 aapp = t*dsqrt( aapp )*d( p )
655 IF( ir1.EQ.0 )notrot = notrot + 1
656 pskipped = pskipped + 1
660 IF( ir1.EQ.0 )notrot = notrot + 1
661 pskipped = pskipped + 1
664 IF( ( i.LE.swband ) .AND.
665 $ ( pskipped.GT.rowskip ) )
THEN
666 IF( ir1.EQ.0 )aapp = -aapp
681 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
682 $ notrot = notrot + min( igl+kbl-1, n ) - p
694 igl = ( ibr-1 )*kbl + 1
696 DO 2010 jbc = ibr + 1, nbl
698 jgl = ( jbc-1 )*kbl + 1
703 DO 2100 p = igl, min( igl+kbl-1, n )
707 IF( aapp.GT.zero )
THEN
711 DO 2200 q = jgl, min( jgl+kbl-1, n )
715 IF( aaqq.GT.zero )
THEN
722 IF( aaqq.GE.one )
THEN
723 IF( aapp.GE.aaqq )
THEN
724 rotok = ( small*aapp ).LE.aaqq
726 rotok = ( small*aaqq ).LE.aapp
728 IF( aapp.LT.( big / aaqq ) )
THEN
729 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
730 $ q ), 1 )*d( p )*d( q ) / aaqq )
733 CALL dcopy( m, a( 1, p ), 1, work, 1 )
734 CALL dlascl(
'G', 0, 0, aapp, d( p ),
735 $ m, 1, work, lda, ierr )
736 aapq =
ddot( m, work, 1, a( 1, q ),
740 IF( aapp.GE.aaqq )
THEN
741 rotok = aapp.LE.( aaqq / small )
743 rotok = aaqq.LE.( aapp / small )
745 IF( aapp.GT.( small / aaqq ) )
THEN
746 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
747 $ q ), 1 )*d( p )*d( q ) / aaqq )
750 CALL dcopy( m, a( 1, q ), 1, work, 1 )
751 CALL dlascl(
'G', 0, 0, aaqq, d( q ),
752 $ m, 1, work, lda, ierr )
753 aapq =
ddot( m, work, 1, a( 1, p ),
758 mxaapq = max( mxaapq, dabs( aapq ) )
762 IF( dabs( aapq ).GT.tol )
THEN
772 theta = -half*dabs( aqoap-apoaq )/aapq
773 IF( aaqq.GT.aapp0 )theta = -theta
775 IF( dabs( theta ).GT.bigtheta )
THEN
777 fastr( 3 ) = t*d( p ) / d( q )
778 fastr( 4 ) = -t*d( q ) / d( p )
779 CALL drotm( m, a( 1, p ), 1,
780 $ a( 1, q ), 1, fastr )
781 IF( rsvec )
CALL drotm( mvl,
785 sva( q ) = aaqq*dsqrt( max( zero,
786 $ one+t*apoaq*aapq ) )
787 aapp = aapp*dsqrt( max( zero,
788 $ one-t*aqoap*aapq ) )
789 mxsinj = max( mxsinj, dabs( t ) )
794 thsign = -dsign( one, aapq )
795 IF( aaqq.GT.aapp0 )thsign = -thsign
796 t = one / ( theta+thsign*
797 $ dsqrt( one+theta*theta ) )
798 cs = dsqrt( one / ( one+t*t ) )
800 mxsinj = max( mxsinj, dabs( sn ) )
801 sva( q ) = aaqq*dsqrt( max( zero,
802 $ one+t*apoaq*aapq ) )
803 aapp = aapp*dsqrt( max( zero,
804 $ one-t*aqoap*aapq ) )
806 apoaq = d( p ) / d( q )
807 aqoap = d( q ) / d( p )
808 IF( d( p ).GE.one )
THEN
810 IF( d( q ).GE.one )
THEN
812 fastr( 4 ) = -t*aqoap
815 CALL drotm( m, a( 1, p ), 1,
818 IF( rsvec )
CALL drotm( mvl,
819 $ v( 1, p ), 1, v( 1, q ),
822 CALL daxpy( m, -t*aqoap,
825 CALL daxpy( m, cs*sn*apoaq,
829 CALL daxpy( mvl, -t*aqoap,
841 IF( d( q ).GE.one )
THEN
842 CALL daxpy( m, t*apoaq,
845 CALL daxpy( m, -cs*sn*aqoap,
849 CALL daxpy( mvl, t*apoaq,
860 IF( d( p ).GE.d( q ) )
THEN
861 CALL daxpy( m, -t*aqoap,
864 CALL daxpy( m, cs*sn*apoaq,
880 CALL daxpy( m, t*apoaq,
891 $ t*apoaq, v( 1, p ),
904 IF( aapp.GT.aaqq )
THEN
905 CALL dcopy( m, a( 1, p ), 1, work,
907 CALL dlascl(
'G', 0, 0, aapp, one,
908 $ m, 1, work, lda, ierr )
909 CALL dlascl(
'G', 0, 0, aaqq, one,
910 $ m, 1, a( 1, q ), lda,
912 temp1 = -aapq*d( p ) / d( q )
913 CALL daxpy( m, temp1, work, 1,
915 CALL dlascl(
'G', 0, 0, one, aaqq,
916 $ m, 1, a( 1, q ), lda,
918 sva( q ) = aaqq*dsqrt( max( zero,
920 mxsinj = max( mxsinj, sfmin )
922 CALL dcopy( m, a( 1, q ), 1, work,
924 CALL dlascl(
'G', 0, 0, aaqq, one,
925 $ m, 1, work, lda, ierr )
926 CALL dlascl(
'G', 0, 0, aapp, one,
927 $ m, 1, a( 1, p ), lda,
929 temp1 = -aapq*d( q ) / d( p )
930 CALL daxpy( m, temp1, work, 1,
932 CALL dlascl(
'G', 0, 0, one, aapp,
933 $ m, 1, a( 1, p ), lda,
935 sva( p ) = aapp*dsqrt( max( zero,
937 mxsinj = max( mxsinj, sfmin )
944 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
946 IF( ( aaqq.LT.rootbig ) .AND.
947 $ ( aaqq.GT.rootsfmin ) )
THEN
948 sva( q ) =
dnrm2( m, a( 1, q ), 1 )*
953 CALL dlassq( m, a( 1, q ), 1, t,
955 sva( q ) = t*dsqrt( aaqq )*d( q )
958 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
959 IF( ( aapp.LT.rootbig ) .AND.
960 $ ( aapp.GT.rootsfmin ) )
THEN
961 aapp =
dnrm2( m, a( 1, p ), 1 )*
966 CALL dlassq( m, a( 1, p ), 1, t,
968 aapp = t*dsqrt( aapp )*d( p )
975 pskipped = pskipped + 1
980 pskipped = pskipped + 1
984 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
990 IF( ( i.LE.swband ) .AND.
991 $ ( pskipped.GT.rowskip ) )
THEN
1004 IF( aapp.EQ.zero )notrot = notrot +
1005 $ min( jgl+kbl-1, n ) - jgl + 1
1006 IF( aapp.LT.zero )notrot = 0
1015 DO 2012 p = igl, min( igl+kbl-1, n )
1016 sva( p ) = dabs( sva( p ) )
1023 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1025 sva( n ) =
dnrm2( m, a( 1, n ), 1 )*d( n )
1029 CALL dlassq( m, a( 1, n ), 1, t, aapp )
1030 sva( n ) = t*dsqrt( aapp )*d( n )
1035 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1036 $ ( iswrot.LE.n ) ) )swband = i
1038 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
1039 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) )
THEN
1043 IF( notrot.GE.emptsw )
GO TO 1994
1060 DO 5991 p = 1, n - 1
1061 q =
idamax( n-p+1, sva( p ), 1 ) + p - 1
1069 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1070 IF( rsvec )
CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
integer function idamax(N, DX, INCX)
IDAMAX
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
subroutine dlassq(N, X, INCX, SCALE, SUMSQ)
DLASSQ updates a sum of squares represented in scaled form.
subroutine drotm(N, DX, INCX, DY, INCY, DPARAM)
DROTM
double precision function dnrm2(N, X, INCX)
DNRM2
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME