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, max, float, min, 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 = min( 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, min( igl+kbl-1, n1 )
406 IF( aapp.GT.zero )
THEN
410 DO 2200 q = jgl, min( 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 = max( 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( max( zero,
485 $ one+t*apoaq*aapq ) )
486 aapp = aapp*sqrt( max( zero,
487 $ one-t*aqoap*aapq ) )
488 mxsinj = max( 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 = max( mxsinj, abs( sn ) )
500 sva( q ) = aaqq*sqrt( max( zero,
501 $ one+t*apoaq*aapq ) )
502 aapp = aapp*sqrt( max( 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( max( zero,
619 mxsinj = max( 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( max( zero,
636 mxsinj = max( 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 $ min( jgl+kbl-1, n ) - jgl + 1
708 IF( aapp.LT.zero )notrot = 0
718 DO 2012 p = igl, min( 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 srotm(N, SX, INCX, SY, INCY, SPARAM)
SROTM
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 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...
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.