236 SUBROUTINE sgsvj1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
237 $ eps, sfmin, tol, nsweep, work, lwork, info )
246 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
250 REAL A( lda, * ), D( n ), SVA( n ), V( ldv, * ),
258 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
261 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
262 $ bigtheta, cs, large, mxaapq, mxsinj, rootbig,
263 $ rooteps, rootsfmin, roottol, small, sn, t,
264 $ temp1, theta, thsign
265 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
266 $ iswrot, jbc, jgl, kbl, mvl, notrot, nblc, nblr,
267 $ p, pskipped, q, rowskip, swband
268 LOGICAL APPLV, ROTOK, RSVEC
274 INTRINSIC abs, amax1, float, min0, sign, sqrt
280 EXTERNAL isamax, lsame, sdot, snrm2
289 applv = lsame( jobv,
'A' )
290 rsvec = lsame( jobv,
'V' )
291 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv,
'N' ) ) )
THEN
293 ELSE IF( m.LT.0 )
THEN
295 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) )
THEN
297 ELSE IF( n1.LT.0 )
THEN
299 ELSE IF( lda.LT.m )
THEN
301 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) )
THEN
303 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304 $ ( applv.AND.( ldv.LT.mv ) ) )
THEN
306 ELSE IF( tol.LE.eps )
THEN
308 ELSE IF( nsweep.LT.0 )
THEN
310 ELSE IF( lwork.LT.m )
THEN
318 CALL
xerbla(
'SGSVJ1', -info )
324 ELSE IF( applv )
THEN
327 rsvec = rsvec .OR. applv
329 rooteps = sqrt( eps )
330 rootsfmin = sqrt( sfmin )
333 rootbig = one / rootsfmin
334 large = big / sqrt( float( m*n ) )
335 bigtheta = one / rooteps
336 roottol = sqrt( tol )
350 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
354 nblc = ( n-n1 ) / kbl
355 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
356 blskip = ( kbl**2 ) + 1
359 rowskip = min0( 5, kbl )
375 DO 1993 i = 1, nsweep
385 DO 2000 ibr = 1, nblr
387 igl = ( ibr-1 )*kbl + 1
393 igl = ( ibr-1 )*kbl + 1
395 DO 2010 jbc = 1, nblc
397 jgl = n1 + ( jbc-1 )*kbl + 1
402 DO 2100 p = igl, min0( igl+kbl-1, n1 )
406 IF( aapp.GT.zero )
THEN
410 DO 2200 q = jgl, min0( jgl+kbl-1, n )
414 IF( aaqq.GT.zero )
THEN
421 IF( aaqq.GE.one )
THEN
422 IF( aapp.GE.aaqq )
THEN
423 rotok = ( small*aapp ).LE.aaqq
425 rotok = ( small*aaqq ).LE.aapp
427 IF( aapp.LT.( big / aaqq ) )
THEN
428 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
429 $ q ), 1 )*d( p )*d( q ) / aaqq )
432 CALL
scopy( m, a( 1, p ), 1, work, 1 )
433 CALL
slascl(
'G', 0, 0, aapp, d( p ),
434 $ m, 1, work, lda, ierr )
435 aapq = sdot( m, work, 1, a( 1, q ),
439 IF( aapp.GE.aaqq )
THEN
440 rotok = aapp.LE.( aaqq / small )
442 rotok = aaqq.LE.( aapp / small )
444 IF( aapp.GT.( small / aaqq ) )
THEN
445 aapq = ( sdot( m, a( 1, p ), 1, a( 1,
446 $ q ), 1 )*d( p )*d( q ) / aaqq )
449 CALL
scopy( m, a( 1, q ), 1, work, 1 )
450 CALL
slascl(
'G', 0, 0, aaqq, d( q ),
451 $ m, 1, work, lda, ierr )
452 aapq = sdot( m, work, 1, a( 1, p ),
457 mxaapq = amax1( mxaapq, abs( aapq ) )
461 IF( abs( aapq ).GT.tol )
THEN
471 theta = -half*abs( aqoap-apoaq ) / aapq
472 IF( aaqq.GT.aapp0 )theta = -theta
474 IF( abs( theta ).GT.bigtheta )
THEN
476 fastr( 3 ) = t*d( p ) / d( q )
477 fastr( 4 ) = -t*d( q ) / d( p )
478 CALL
srotm( m, a( 1, p ), 1,
479 $ a( 1, q ), 1, fastr )
480 IF( rsvec )CALL
srotm( mvl,
484 sva( q ) = aaqq*sqrt( amax1( zero,
485 $ one+t*apoaq*aapq ) )
486 aapp = aapp*sqrt( amax1( zero,
487 $ one-t*aqoap*aapq ) )
488 mxsinj = amax1( mxsinj, abs( t ) )
493 thsign = -sign( one, aapq )
494 IF( aaqq.GT.aapp0 )thsign = -thsign
495 t = one / ( theta+thsign*
496 $ sqrt( one+theta*theta ) )
497 cs = sqrt( one / ( one+t*t ) )
499 mxsinj = amax1( mxsinj, abs( sn ) )
500 sva( q ) = aaqq*sqrt( amax1( zero,
501 $ one+t*apoaq*aapq ) )
502 aapp = aapp*sqrt( amax1( zero,
503 $ one-t*aqoap*aapq ) )
505 apoaq = d( p ) / d( q )
506 aqoap = d( q ) / d( p )
507 IF( d( p ).GE.one )
THEN
509 IF( d( q ).GE.one )
THEN
511 fastr( 4 ) = -t*aqoap
514 CALL
srotm( m, a( 1, p ), 1,
517 IF( rsvec )CALL
srotm( mvl,
518 $ v( 1, p ), 1, v( 1, q ),
521 CALL
saxpy( m, -t*aqoap,
524 CALL
saxpy( m, cs*sn*apoaq,
528 CALL
saxpy( mvl, -t*aqoap,
540 IF( d( q ).GE.one )
THEN
541 CALL
saxpy( m, t*apoaq,
544 CALL
saxpy( m, -cs*sn*aqoap,
548 CALL
saxpy( mvl, t*apoaq,
559 IF( d( p ).GE.d( q ) )
THEN
560 CALL
saxpy( m, -t*aqoap,
563 CALL
saxpy( m, cs*sn*apoaq,
579 CALL
saxpy( m, t*apoaq,
590 $ t*apoaq, v( 1, p ),
603 IF( aapp.GT.aaqq )
THEN
604 CALL
scopy( m, a( 1, p ), 1, work,
606 CALL
slascl(
'G', 0, 0, aapp, one,
607 $ m, 1, work, lda, ierr )
608 CALL
slascl(
'G', 0, 0, aaqq, one,
609 $ m, 1, a( 1, q ), lda,
611 temp1 = -aapq*d( p ) / d( q )
612 CALL
saxpy( m, temp1, work, 1,
614 CALL
slascl(
'G', 0, 0, one, aaqq,
615 $ m, 1, a( 1, q ), lda,
617 sva( q ) = aaqq*sqrt( amax1( zero,
619 mxsinj = amax1( mxsinj, sfmin )
621 CALL
scopy( m, a( 1, q ), 1, work,
623 CALL
slascl(
'G', 0, 0, aaqq, one,
624 $ m, 1, work, lda, ierr )
625 CALL
slascl(
'G', 0, 0, aapp, one,
626 $ m, 1, a( 1, p ), lda,
628 temp1 = -aapq*d( q ) / d( p )
629 CALL
saxpy( m, temp1, work, 1,
631 CALL
slascl(
'G', 0, 0, one, aapp,
632 $ m, 1, a( 1, p ), lda,
634 sva( p ) = aapp*sqrt( amax1( zero,
636 mxsinj = amax1( mxsinj, sfmin )
643 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
645 IF( ( aaqq.LT.rootbig ) .AND.
646 $ ( aaqq.GT.rootsfmin ) )
THEN
647 sva( q ) = snrm2( m, a( 1, q ), 1 )*
652 CALL
slassq( m, a( 1, q ), 1, t,
654 sva( q ) = t*sqrt( aaqq )*d( q )
657 IF( ( aapp / aapp0 )**2.LE.rooteps )
THEN
658 IF( ( aapp.LT.rootbig ) .AND.
659 $ ( aapp.GT.rootsfmin ) )
THEN
660 aapp = snrm2( m, a( 1, p ), 1 )*
665 CALL
slassq( m, a( 1, p ), 1, t,
667 aapp = t*sqrt( aapp )*d( p )
675 pskipped = pskipped + 1
680 pskipped = pskipped + 1
685 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
691 IF( ( i.LE.swband ) .AND.
692 $ ( pskipped.GT.rowskip ) )
THEN
706 IF( aapp.EQ.zero )notrot = notrot +
707 $ min0( jgl+kbl-1, n ) - jgl + 1
708 IF( aapp.LT.zero )notrot = 0
718 DO 2012 p = igl, min0( igl+kbl-1, n )
719 sva( p ) = abs( sva( p ) )
726 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
728 sva( n ) = snrm2( m, a( 1, n ), 1 )*d( n )
732 CALL
slassq( m, a( 1, n ), 1, t, aapp )
733 sva( n ) = t*sqrt( aapp )*d( n )
738 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
739 $ ( iswrot.LE.n ) ) )swband = i
741 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.float( n )*tol ) .AND.
742 $ ( float( n )*mxaapq*mxsinj.LT.tol ) )
THEN
747 IF( notrot.GE.emptsw )go to 1994
766 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
774 CALL
sswap( m, a( 1, p ), 1, a( 1, q ), 1 )
775 IF( rsvec )CALL
sswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
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 sgsvj1(JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV, EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO)
SGSVJ1 pre-processor for the routine sgesvj, applies Jacobi rotations targeting only particular pivot...