218 SUBROUTINE sgsvj0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS,
219 $ sfmin, tol, nsweep, work, lwork, info )
227 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
232 REAL A( lda, * ), SVA( n ), D( n ), V( ldv, * ),
240 parameter ( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
243 REAL 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
256 INTRINSIC abs, max, float, min, sign, sqrt
262 EXTERNAL isamax, lsame, sdot, snrm2
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(
'SGSVJ0', -info )
304 ELSE IF( applv )
THEN
307 rsvec = rsvec .OR. applv
309 rooteps = sqrt( eps )
310 rootsfmin = sqrt( sfmin )
313 rootbig = one / rootsfmin
314 bigtheta = one / rooteps
315 roottol = sqrt( 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 = isamax( n-p+1, sva( p ), 1 ) + p - 1
375 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
376 IF( rsvec )
CALL sswap( mvl, v( 1, p ), 1,
400 IF( ( sva( p ).LT.rootbig ) .AND.
401 $ ( sva( p ).GT.rootsfmin ) )
THEN
402 sva( p ) = snrm2( m, a( 1, p ), 1 )*d( p )
406 CALL slassq( m, a( 1, p ), 1, temp1, aapp )
407 sva( p ) = temp1*sqrt( 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 = ( sdot( m, a( 1, p ), 1, a( 1,
430 $ q ), 1 )*d( p )*d( q ) / aaqq )
433 CALL scopy( m, a( 1, p ), 1, work, 1 )
434 CALL slascl(
'G', 0, 0, aapp, d( p ),
435 $ m, 1, work, lda, ierr )
436 aapq = sdot( m, work, 1, a( 1, q ),
440 rotok = aapp.LE.( aaqq / small )
441 IF( aapp.GT.( small / aaqq ) )
THEN
442 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
443 $ q ), 1 )*d( p )*d( q ) / aaqq )
446 CALL scopy( m, a( 1, q ), 1, work, 1 )
447 CALL slascl(
'G', 0, 0, aaqq, d( q ),
448 $ m, 1, work, lda, ierr )
449 aapq = sdot( m, work, 1, a( 1, p ),
454 mxaapq = max( mxaapq, abs( aapq ) )
458 IF( abs( aapq ).GT.tol )
THEN
473 theta = -half*abs( aqoap-apoaq ) / aapq
475 IF( abs( theta ).GT.bigtheta )
THEN
478 fastr( 3 ) = t*d( p ) / d( q )
479 fastr( 4 ) = -t*d( q ) / d( p )
480 CALL srotm( m, a( 1, p ), 1,
481 $ a( 1, q ), 1, fastr )
482 IF( rsvec )
CALL srotm( mvl,
486 sva( q ) = aaqq*sqrt( max( zero,
487 $ one+t*apoaq*aapq ) )
488 aapp = aapp*sqrt( max( zero,
489 $ one-t*aqoap*aapq ) )
490 mxsinj = max( mxsinj, abs( t ) )
496 thsign = -sign( one, aapq )
497 t = one / ( theta+thsign*
498 $ sqrt( one+theta*theta ) )
499 cs = sqrt( one / ( one+t*t ) )
502 mxsinj = max( mxsinj, abs( sn ) )
503 sva( q ) = aaqq*sqrt( max( zero,
504 $ one+t*apoaq*aapq ) )
505 aapp = aapp*sqrt( 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 srotm( m, a( 1, p ), 1,
519 IF( rsvec )
CALL srotm( mvl,
520 $ v( 1, p ), 1, v( 1, q ),
523 CALL saxpy( m, -t*aqoap,
526 CALL saxpy( m, cs*sn*apoaq,
532 CALL saxpy( mvl, -t*aqoap,
542 IF( d( q ).GE.one )
THEN
543 CALL saxpy( m, t*apoaq,
546 CALL saxpy( m, -cs*sn*aqoap,
552 CALL saxpy( mvl, t*apoaq,
561 IF( d( p ).GE.d( q ) )
THEN
562 CALL saxpy( m, -t*aqoap,
565 CALL saxpy( m, cs*sn*apoaq,
581 CALL saxpy( m, t*apoaq,
592 $ t*apoaq, v( 1, p ),
606 CALL scopy( m, a( 1, p ), 1, work, 1 )
607 CALL slascl(
'G', 0, 0, aapp, one, m,
608 $ 1, work, lda, ierr )
609 CALL slascl(
'G', 0, 0, aaqq, one, m,
610 $ 1, a( 1, q ), lda, ierr )
611 temp1 = -aapq*d( p ) / d( q )
612 CALL saxpy( m, temp1, work, 1,
614 CALL slascl(
'G', 0, 0, one, aaqq, m,
615 $ 1, a( 1, q ), lda, ierr )
616 sva( q ) = aaqq*sqrt( 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 ) = snrm2( m, a( 1, q ), 1 )*
633 CALL slassq( m, a( 1, q ), 1, t,
635 sva( q ) = t*sqrt( aaqq )*d( q )
638 IF( ( aapp / aapp0 ).LE.rooteps )
THEN
639 IF( ( aapp.LT.rootbig ) .AND.
640 $ ( aapp.GT.rootsfmin ) )
THEN
641 aapp = snrm2( m, a( 1, p ), 1 )*
646 CALL slassq( m, a( 1, p ), 1, t,
648 aapp = t*sqrt( 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 = ( sdot( m, a( 1, p ), 1, a( 1,
730 $ q ), 1 )*d( p )*d( q ) / aaqq )
733 CALL scopy( m, a( 1, p ), 1, work, 1 )
734 CALL slascl(
'G', 0, 0, aapp, d( p ),
735 $ m, 1, work, lda, ierr )
736 aapq = sdot( 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 = ( sdot( m, a( 1, p ), 1, a( 1,
747 $ q ), 1 )*d( p )*d( q ) / aaqq )
750 CALL scopy( m, a( 1, q ), 1, work, 1 )
751 CALL slascl(
'G', 0, 0, aaqq, d( q ),
752 $ m, 1, work, lda, ierr )
753 aapq = sdot( m, work, 1, a( 1, p ),
758 mxaapq = max( mxaapq, abs( aapq ) )
762 IF( abs( aapq ).GT.tol )
THEN
772 theta = -half*abs( aqoap-apoaq ) / aapq
773 IF( aaqq.GT.aapp0 )theta = -theta
775 IF( abs( theta ).GT.bigtheta )
THEN
777 fastr( 3 ) = t*d( p ) / d( q )
778 fastr( 4 ) = -t*d( q ) / d( p )
779 CALL srotm( m, a( 1, p ), 1,
780 $ a( 1, q ), 1, fastr )
781 IF( rsvec )
CALL srotm( mvl,
785 sva( q ) = aaqq*sqrt( max( zero,
786 $ one+t*apoaq*aapq ) )
787 aapp = aapp*sqrt( max( zero,
788 $ one-t*aqoap*aapq ) )
789 mxsinj = max( mxsinj, abs( t ) )
794 thsign = -sign( one, aapq )
795 IF( aaqq.GT.aapp0 )thsign = -thsign
796 t = one / ( theta+thsign*
797 $ sqrt( one+theta*theta ) )
798 cs = sqrt( one / ( one+t*t ) )
800 mxsinj = max( mxsinj, abs( sn ) )
801 sva( q ) = aaqq*sqrt( max( zero,
802 $ one+t*apoaq*aapq ) )
803 aapp = aapp*sqrt( 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 srotm( m, a( 1, p ), 1,
818 IF( rsvec )
CALL srotm( mvl,
819 $ v( 1, p ), 1, v( 1, q ),
822 CALL saxpy( m, -t*aqoap,
825 CALL saxpy( m, cs*sn*apoaq,
829 CALL saxpy( mvl, -t*aqoap,
841 IF( d( q ).GE.one )
THEN
842 CALL saxpy( m, t*apoaq,
845 CALL saxpy( m, -cs*sn*aqoap,
849 CALL saxpy( mvl, t*apoaq,
860 IF( d( p ).GE.d( q ) )
THEN
861 CALL saxpy( m, -t*aqoap,
864 CALL saxpy( m, cs*sn*apoaq,
880 CALL saxpy( m, t*apoaq,
891 $ t*apoaq, v( 1, p ),
904 IF( aapp.GT.aaqq )
THEN
905 CALL scopy( m, a( 1, p ), 1, work,
907 CALL slascl(
'G', 0, 0, aapp, one,
908 $ m, 1, work, lda, ierr )
909 CALL slascl(
'G', 0, 0, aaqq, one,
910 $ m, 1, a( 1, q ), lda,
912 temp1 = -aapq*d( p ) / d( q )
913 CALL saxpy( m, temp1, work, 1,
915 CALL slascl(
'G', 0, 0, one, aaqq,
916 $ m, 1, a( 1, q ), lda,
918 sva( q ) = aaqq*sqrt( max( zero,
920 mxsinj = max( mxsinj, sfmin )
922 CALL scopy( m, a( 1, q ), 1, work,
924 CALL slascl(
'G', 0, 0, aaqq, one,
925 $ m, 1, work, lda, ierr )
926 CALL slascl(
'G', 0, 0, aapp, one,
927 $ m, 1, a( 1, p ), lda,
929 temp1 = -aapq*d( q ) / d( p )
930 CALL saxpy( m, temp1, work, 1,
932 CALL slascl(
'G', 0, 0, one, aapp,
933 $ m, 1, a( 1, p ), lda,
935 sva( p ) = aapp*sqrt( 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 ) = snrm2( m, a( 1, q ), 1 )*
953 CALL slassq( m, a( 1, q ), 1, t,
955 sva( q ) = t*sqrt( aaqq )*d( q )
958 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
959 IF( ( aapp.LT.rootbig ) .AND.
960 $ ( aapp.GT.rootsfmin ) )
THEN
961 aapp = snrm2( m, a( 1, p ), 1 )*
966 CALL slassq( m, a( 1, p ), 1, t,
968 aapp = t*sqrt( 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 ) = abs( sva( p ) )
1023 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1025 sva( n ) = snrm2( m, a( 1, n ), 1 )*d( n )
1029 CALL slassq( m, a( 1, n ), 1, t, aapp )
1030 sva( n ) = t*sqrt( 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.float( n )*tol ) .AND.
1039 $ ( float( n )*mxaapq*mxsinj.LT.tol ) )
THEN
1043 IF( notrot.GE.emptsw )
GO TO 1994
1060 DO 5991 p = 1, n - 1
1061 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
1069 CALL sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1070 IF( rsvec )
CALL sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
subroutine sgsvj0(JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
SGSVJ0 pre-processor for the routine sgesvj.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.