226 SUBROUTINE sbdsvdx( UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU,
227 $ ns, s, z, ldz, work, iwork, info)
235 CHARACTER JOBZ, RANGE, UPLO
236 INTEGER IL, INFO, IU, LDZ, N, NS
241 REAL D( * ), E( * ), S( * ), WORK( * ),
248 REAL ZERO, ONE, TEN, HNDRD, MEIGTH
249 parameter ( zero = 0.0e0, one = 1.0e0, ten = 10.0e0,
250 $ hndrd = 100.0e0, meigth = -0.1250e0 )
252 parameter ( fudge = 2.0e0 )
256 LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
257 INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
258 $ ietgk, iifail, iiwork, iltgk, irowu, irowv,
259 $ irowz, isbeg, isplt, itemp, iutgk, j, k,
260 $ ntgk, nru, nrv, nsl
261 REAL ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
262 $ smin, sqrt2, thresh, tol, ulp,
263 $ vltgk, vutgk, zjtji
268 REAL SDOT, SLAMCH, SNRM2
269 EXTERNAL isamax, lsame,
saxpy, sdot, slamch, snrm2
275 INTRINSIC abs,
REAL, SIGN, SQRT
281 allsv = lsame( range,
'A' )
282 valsv = lsame( range,
'V' )
283 indsv = lsame( range,
'I' )
284 wantz = lsame( jobz,
'V' )
285 lower = lsame( uplo,
'L' )
288 IF( .NOT.lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
290 ELSE IF( .NOT.( wantz .OR. lsame( jobz,
'N' ) ) )
THEN
292 ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) )
THEN
294 ELSE IF( n.LT.0 )
THEN
296 ELSE IF( n.GT.0 )
THEN
298 IF( vl.LT.zero )
THEN
300 ELSE IF( vu.LE.vl )
THEN
303 ELSE IF( indsv )
THEN
304 IF( il.LT.1 .OR. il.GT.max( 1, n ) )
THEN
306 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n )
THEN
312 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
316 CALL xerbla(
'SBDSVDX', -info )
326 IF( allsv .OR. indsv )
THEN
328 s( 1 ) = abs( d( 1 ) )
330 IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) )
THEN
332 s( 1 ) = abs( d( 1 ) )
336 z( 1, 1 ) = sign( one, d( 1 ) )
342 abstol = 2*slamch(
'Safe Minimum' )
343 ulp = slamch(
'Precision' )
344 eps = slamch(
'Epsilon' )
345 sqrt2 = sqrt( 2.0e0 )
353 tol = max( ten, min( hndrd, eps**meigth ) )*eps
357 i = isamax( n, d, 1 )
359 i = isamax( n-1, e, 1 )
360 smax = max( smax, abs( e( i ) ) )
365 IF( smin.NE.zero )
THEN
368 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
369 smin = min( smin, mu )
370 IF( smin.EQ.zero )
EXIT
373 smin = smin / sqrt(
REAL( N ) )
379 IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
380 IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
382 IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
390 iiwork = iifail + n*2
409 IF( wantz )
CALL slaset(
'F', n*2, n+1, zero, zero, z, ldz )
410 ELSE IF( valsv )
THEN
419 work( idtgk:idtgk+2*n-1 ) = zero
420 CALL scopy( n, d, 1, work( ietgk ), 2 )
421 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
422 CALL sstevx(
'N',
'V', n*2, work( idtgk ), work( ietgk ),
423 $ vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
424 $ z, ldz, work( itemp ), iwork( iiwork ),
425 $ iwork( iifail ), info )
429 IF( wantz )
CALL slaset(
'F', n*2, ns, zero, zero, z, ldz )
431 ELSE IF( indsv )
THEN
443 work( idtgk:idtgk+2*n-1 ) = zero
444 CALL scopy( n, d, 1, work( ietgk ), 2 )
445 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
446 CALL sstevx(
'N',
'I', n*2, work( idtgk ), work( ietgk ),
447 $ vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
448 $ z, ldz, work( itemp ), iwork( iiwork ),
449 $ iwork( iifail ), info )
450 vltgk = s( 1 ) - fudge*smax*ulp*n
451 work( idtgk:idtgk+2*n-1 ) = zero
452 CALL scopy( n, d, 1, work( ietgk ), 2 )
453 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
454 CALL sstevx(
'N',
'I', n*2, work( idtgk ), work( ietgk ),
455 $ vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
456 $ z, ldz, work( itemp ), iwork( iiwork ),
457 $ iwork( iifail ), info )
458 vutgk = s( 1 ) + fudge*smax*ulp*n
459 vutgk = min( vutgk, zero )
464 IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
466 IF( wantz )
CALL slaset(
'F', n*2, iu-il+1, zero, zero, z, ldz)
491 work( ietgk+2*n-1 ) = zero
492 work( idtgk:idtgk+2*n-1 ) = zero
493 CALL scopy( n, d, 1, work( ietgk ), 2 )
494 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
501 IF( work( ietgk+ieptr-1 ).EQ.zero )
THEN
508 DO idptr = idbeg, idend, 2
509 IF( work( ietgk+idptr-1 ).EQ.zero )
THEN
514 IF( idptr.EQ.idbeg )
THEN
519 IF( idbeg.EQ.idend)
THEN
523 ELSE IF( idptr.EQ.idend )
THEN
528 nru = (idend-isplt)/2 + 1
530 IF( isplt.NE.idbeg )
THEN
534 IF( isplt.EQ.idbeg )
THEN
538 nru = (idptr-idbeg)/2
544 nru = (idptr-isplt)/2 + 1
548 ELSE IF( idptr.EQ.idend )
THEN
552 IF( isplt.EQ.idbeg )
THEN
556 nru = (idend-idbeg)/2 + 1
562 nrv = (idend-isplt)/2 + 1
579 IF( allsv .OR. vutgk.EQ.zero )
THEN
582 $ mod(ntgk,2).GT.0 )
THEN
593 CALL sstevx( jobz, rngvx, ntgk, work( idtgk+isplt-1 ),
603 emin = abs( maxval( s( isbeg:isbeg+nsl-1 ) ) )
605 IF( nsl.GT.0 .AND. wantz )
THEN
616 $ vutgk.EQ.zero .AND.
617 $ mod(ntgk,2).EQ.0 .AND.
618 $ emin.EQ.0 .AND. .NOT.split )
THEN
625 z( irowz:irowz+ntgk-1,icolz+nsl-2 ) =
626 $ z( irowz:irowz+ntgk-1,icolz+nsl-2 ) +
627 $ z( irowz:irowz+ntgk-1,icolz+nsl-1 )
628 z( irowz:irowz+ntgk-1,icolz+nsl-1 ) =
636 DO i = 0, min( nsl-1, nru-1 )
637 nrmu = snrm2( nru, z( irowu, icolz+i ), 2 )
638 IF( nrmu.EQ.zero )
THEN
642 CALL sscal( nru, one/nrmu,
643 $ z( irowu,icolz+i ), 2 )
644 IF( nrmu.NE.one .AND.
645 $ abs( nrmu-ortol )*sqrt2.GT.one )
648 zjtji = -sdot( nru, z( irowu, icolz+j ),
649 $ 2, z( irowu, icolz+i ), 2 )
650 CALL saxpy( nru, zjtji,
651 $ z( irowu, icolz+j ), 2,
652 $ z( irowu, icolz+i ), 2 )
654 nrmu = snrm2( nru, z( irowu, icolz+i ), 2 )
655 CALL sscal( nru, one/nrmu,
656 $ z( irowu,icolz+i ), 2 )
659 DO i = 0, min( nsl-1, nrv-1 )
660 nrmv = snrm2( nrv, z( irowv, icolz+i ), 2 )
661 IF( nrmv.EQ.zero )
THEN
665 CALL sscal( nrv, -one/nrmv,
666 $ z( irowv,icolz+i ), 2 )
667 IF( nrmv.NE.one .AND.
668 $ abs( nrmv-ortol )*sqrt2.GT.one )
671 zjtji = -sdot( nrv, z( irowv, icolz+j ),
672 $ 2, z( irowv, icolz+i ), 2 )
673 CALL saxpy( nru, zjtji,
674 $ z( irowv, icolz+j ), 2,
675 $ z( irowv, icolz+i ), 2 )
677 nrmv = snrm2( nrv, z( irowv, icolz+i ), 2 )
678 CALL sscal( nrv, one/nrmv,
679 $ z( irowv,icolz+i ), 2 )
682 IF( vutgk.EQ.zero .AND.
683 $ idptr.LT.idend .AND.
684 $ mod(ntgk,2).GT.0 )
THEN
692 z( irowz:irowz+ntgk-1,n+1 ) =
693 $ z( irowz:irowz+ntgk-1,ns+nsl )
694 z( irowz:irowz+ntgk-1,ns+nsl ) =
699 nsl = min( nsl, nru )
705 s( isbeg+i ) = abs( s( isbeg+i ) )
720 IF( irowz.LT.n*2 .AND. wantz )
THEN
721 z( 1:irowz-1, icolz ) = zero
724 IF( split .AND. wantz )
THEN
729 z( idbeg:idend-ntgk+1,isbeg-1 ) =
730 $ z( idbeg:idend-ntgk+1,isbeg-1 ) +
731 $ z( idbeg:idend-ntgk+1,n+1 )
732 z( idbeg:idend-ntgk+1,n+1 ) = 0
749 IF( s( j ).LE.smin )
THEN
754 IF( k.NE.ns+1-i )
THEN
757 IF( wantz )
CALL sswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ), 1 )
767 IF( wantz ) z( 1:n*2,k+1:ns ) = zero
777 CALL scopy( n*2, z( 1,i ), 1, work, 1 )
779 CALL scopy( n, work( 2 ), 2, z( n+1,i ), 1 )
780 CALL scopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
782 CALL scopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
783 CALL scopy( n, work( 1 ), 2, z( n+1,i ), 1 )
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
subroutine sstevx(JOBZ, RANGE, N, D, E, VL, VU, IL, IU, ABSTOL, M, W, Z, LDZ, WORK, IWORK, IFAIL, INFO)
SSTEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for OTHER matric...
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine sswap(N, SX, INCX, SY, INCY)
SSWAP
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY