226
227
228
229
230
231
232 CHARACTER JOBZ, RANGE, UPLO
233 INTEGER IL, INFO, IU, LDZ, N, NS
234 REAL VL, VU
235
236
237 INTEGER IWORK( * )
238 REAL D( * ), E( * ), S( * ), WORK( * ),
239 $ Z( LDZ, * )
240
241
242
243
244
245 REAL ZERO, ONE, TEN, HNDRD, MEIGTH
246 parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0,
247 $ hndrd = 100.0e0, meigth = -0.1250e0 )
248 REAL FUDGE
249 parameter( fudge = 2.0e0 )
250
251
252 CHARACTER RNGVX
253 LOGICAL ALLSV, INDSV, LOWER, SPLIT, SVEQ0, VALSV, WANTZ
254 INTEGER I, ICOLZ, IDBEG, IDEND, IDTGK, IDPTR, IEPTR,
255 $ IETGK, IIFAIL, IIWORK, ILTGK, IROWU, IROWV,
256 $ IROWZ, ISBEG, ISPLT, ITEMP, IUTGK, J, K,
257 $ NTGK, NRU, NRV, NSL
258 REAL ABSTOL, EPS, EMIN, MU, NRMU, NRMV, ORTOL, SMAX,
259 $ SMIN, SQRT2, THRESH, TOL, ULP,
260 $ VLTGK, VUTGK, ZJTJI
261
262
263 LOGICAL LSAME
264 INTEGER ISAMAX
265 REAL SDOT, SLAMCH, SNRM2
267
268
270
271
272 INTRINSIC abs, real, sign, sqrt
273
274
275
276
277
278 allsv =
lsame( range,
'A' )
279 valsv =
lsame( range,
'V' )
280 indsv =
lsame( range,
'I' )
281 wantz =
lsame( jobz,
'V' )
282 lower =
lsame( uplo,
'L' )
283
284 info = 0
285 IF( .NOT.
lsame( uplo,
'U' ) .AND. .NOT.lower )
THEN
286 info = -1
287 ELSE IF( .NOT.( wantz .OR.
lsame( jobz,
'N' ) ) )
THEN
288 info = -2
289 ELSE IF( .NOT.( allsv .OR. valsv .OR. indsv ) ) THEN
290 info = -3
291 ELSE IF( n.LT.0 ) THEN
292 info = -4
293 ELSE IF( n.GT.0 ) THEN
294 IF( valsv ) THEN
295 IF( vl.LT.zero ) THEN
296 info = -7
297 ELSE IF( vu.LE.vl ) THEN
298 info = -8
299 END IF
300 ELSE IF( indsv ) THEN
301 IF( il.LT.1 .OR. il.GT.max( 1, n ) ) THEN
302 info = -9
303 ELSE IF( iu.LT.min( n, il ) .OR. iu.GT.n ) THEN
304 info = -10
305 END IF
306 END IF
307 END IF
308 IF( info.EQ.0 ) THEN
309 IF( ldz.LT.1 .OR. ( wantz .AND. ldz.LT.n*2 ) ) info = -14
310 END IF
311
312 IF( info.NE.0 ) THEN
313 CALL xerbla(
'SBDSVDX', -info )
314 RETURN
315 END IF
316
317
318
319 ns = 0
320 IF( n.EQ.0 ) RETURN
321
322 IF( n.EQ.1 ) THEN
323 IF( allsv .OR. indsv ) THEN
324 ns = 1
325 s( 1 ) = abs( d( 1 ) )
326 ELSE
327 IF( vl.LT.abs( d( 1 ) ) .AND. vu.GE.abs( d( 1 ) ) ) THEN
328 ns = 1
329 s( 1 ) = abs( d( 1 ) )
330 END IF
331 END IF
332 IF( wantz ) THEN
333 z( 1, 1 ) = sign( one, d( 1 ) )
334 z( 2, 1 ) = one
335 END IF
336 RETURN
337 END IF
338
339 abstol = 2*
slamch(
'Safe Minimum' )
340 ulp =
slamch(
'Precision' )
342 sqrt2 = sqrt( 2.0e0 )
343 ortol = sqrt( ulp )
344
345
346
347
348
349
350 tol = max( ten, min( hndrd, eps**meigth ) )*eps
351
352
353
355 smax = abs( d( i ) )
357 smax = max( smax, abs( e( i ) ) )
358
359
360
361 smin = abs( d( 1 ) )
362 IF( smin.NE.zero ) THEN
363 mu = smin
364 DO i = 2, n
365 mu = abs( d( i ) )*( mu / ( mu+abs( e( i-1 ) ) ) )
366 smin = min( smin, mu )
367 IF( smin.EQ.zero ) EXIT
368 END DO
369 END IF
370 smin = smin / sqrt( real( n ) )
371 thresh = tol*smin
372
373
374
375 DO i = 1, n-1
376 IF( abs( d( i ) ).LE.thresh ) d( i ) = zero
377 IF( abs( e( i ) ).LE.thresh ) e( i ) = zero
378 END DO
379 IF( abs( d( n ) ).LE.thresh ) d( n ) = zero
380
381
382
383 idtgk = 1
384 ietgk = idtgk + n*2
385 itemp = ietgk + n*2
386 iifail = 1
387 iiwork = iifail + n*2
388
389
390
391
392
393 iltgk = 0
394 iutgk = 0
395 vltgk = zero
396 vutgk = zero
397
398 IF( allsv ) THEN
399
400
401
402
403
404
405 rngvx = 'I'
406 IF( wantz )
CALL slaset(
'F', n*2, n+1, zero, zero, z, ldz )
407 ELSE IF( valsv ) THEN
408
409
410
411
412
413 rngvx = 'V'
414 vltgk = -vu
415 vutgk = -vl
416 work( idtgk:idtgk+2*n-1 ) = zero
417 CALL scopy( n, d, 1, work( ietgk ), 2 )
418 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
419 CALL sstevx(
'N',
'V', n*2, work( idtgk ), work( ietgk ),
420 $ vltgk, vutgk, iltgk, iltgk, abstol, ns, s,
421 $ z, ldz, work( itemp ), iwork( iiwork ),
422 $ iwork( iifail ), info )
423 IF( ns.EQ.0 ) THEN
424 RETURN
425 ELSE
426 IF( wantz )
CALL slaset(
'F', n*2, ns, zero, zero, z, ldz )
427 END IF
428 ELSE IF( indsv ) THEN
429
430
431
432
433
434
435
436
437 iltgk = il
438 iutgk = iu
439 rngvx = 'V'
440 work( idtgk:idtgk+2*n-1 ) = zero
441 CALL scopy( n, d, 1, work( ietgk ), 2 )
442 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
443 CALL sstevx(
'N',
'I', n*2, work( idtgk ), work( ietgk ),
444 $ vltgk, vltgk, iltgk, iltgk, abstol, ns, s,
445 $ z, ldz, work( itemp ), iwork( iiwork ),
446 $ iwork( iifail ), info )
447 vltgk = s( 1 ) - fudge*smax*ulp*n
448 work( idtgk:idtgk+2*n-1 ) = zero
449 CALL scopy( n, d, 1, work( ietgk ), 2 )
450 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
451 CALL sstevx(
'N',
'I', n*2, work( idtgk ), work( ietgk ),
452 $ vutgk, vutgk, iutgk, iutgk, abstol, ns, s,
453 $ z, ldz, work( itemp ), iwork( iiwork ),
454 $ iwork( iifail ), info )
455 vutgk = s( 1 ) + fudge*smax*ulp*n
456 vutgk = min( vutgk, zero )
457
458
459
460
461 IF( vltgk.EQ.vutgk ) vltgk = vltgk - tol
462
463 IF( wantz )
CALL slaset(
'F', n*2, iu-il+1, zero, zero, z, ldz)
464 END IF
465
466
467
468
469
470
471
472
473 ns = 0
474 nru = 0
475 nrv = 0
476 idbeg = 1
477 isbeg = 1
478 irowz = 1
479 icolz = 1
480 irowu = 2
481 irowv = 1
482 split = .false.
483 sveq0 = .false.
484
485
486
487 s( 1:n ) = zero
488 work( ietgk+2*n-1 ) = zero
489 work( idtgk:idtgk+2*n-1 ) = zero
490 CALL scopy( n, d, 1, work( ietgk ), 2 )
491 CALL scopy( n-1, e, 1, work( ietgk+1 ), 2 )
492
493
494
495
496
497 DO ieptr = 2, n*2, 2
498 IF( work( ietgk+ieptr-1 ).EQ.zero ) THEN
499
500
501
502
503 isplt = idbeg
504 idend = ieptr - 1
505 DO idptr = idbeg, idend, 2
506 IF( work( ietgk+idptr-1 ).EQ.zero ) THEN
507
508
509
510
511 IF( idptr.EQ.idbeg ) THEN
512
513
514
515 sveq0 = .true.
516 IF( idbeg.EQ.idend) THEN
517 nru = 1
518 nrv = 1
519 END IF
520 ELSE IF( idptr.EQ.idend ) THEN
521
522
523
524 sveq0 = .true.
525 nru = (idend-isplt)/2 + 1
526 nrv = nru
527 IF( isplt.NE.idbeg ) THEN
528 nru = nru + 1
529 END IF
530 ELSE
531 IF( isplt.EQ.idbeg ) THEN
532
533
534
535 nru = (idptr-idbeg)/2
536 nrv = nru + 1
537 ELSE
538
539
540
541 nru = (idptr-isplt)/2 + 1
542 nrv = nru
543 END IF
544 END IF
545 ELSE IF( idptr.EQ.idend ) THEN
546
547
548
549 IF( isplt.EQ.idbeg ) THEN
550
551
552
553 nru = (idend-idbeg)/2 + 1
554 nrv = nru
555 ELSE
556
557
558
559 nrv = (idend-isplt)/2 + 1
560 nru = nrv + 1
561 END IF
562 END IF
563
564 ntgk = nru + nrv
565
566 IF( ntgk.GT.0 ) THEN
567
568
569
570
571
572
573
574 iltgk = 1
575 iutgk = ntgk / 2
576 IF( allsv .OR. vutgk.EQ.zero ) THEN
577 IF( sveq0 .OR.
578 $ smin.LT.eps .OR.
579 $ mod(ntgk,2).GT.0 ) THEN
580
581
582 iutgk = iutgk + 1
583 END IF
584 END IF
585
586
587
588
589
590 CALL sstevx( jobz, rngvx, ntgk, work( idtgk+isplt-1 ),
591 $ work( ietgk+isplt-1 ), vltgk, vutgk,
592 $ iltgk, iutgk, abstol, nsl, s( isbeg ),
593 $ z( irowz,icolz ), ldz, work( itemp ),
594 $ iwork( iiwork ), iwork( iifail ),
595 $ info )
596 IF( info.NE.0 ) THEN
597
598 RETURN
599 END IF
600 emin = abs( maxval( s( isbeg:isbeg+nsl-1 ) ) )
601
602 IF( nsl.GT.0 .AND. wantz ) THEN
603
604
605
606
607
608
609
610
611
612 IF( nsl.GT.1 .AND.
613 $ vutgk.EQ.zero .AND.
614 $ mod(ntgk,2).EQ.0 .AND.
615 $ emin.EQ.0 .AND. .NOT.split ) THEN
616
617
618
619
620
621
622 z( irowz:irowz+ntgk-1,icolz+nsl-2 ) =
623 $ z( irowz:irowz+ntgk-1,icolz+nsl-2 ) +
624 $ z( irowz:irowz+ntgk-1,icolz+nsl-1 )
625 z( irowz:irowz+ntgk-1,icolz+nsl-1 ) =
626 $ zero
627
628
629
630
631 END IF
632
633 DO i = 0, min( nsl-1, nru-1 )
634 nrmu =
snrm2( nru, z( irowu, icolz+i ), 2 )
635 IF( nrmu.EQ.zero ) THEN
636 info = n*2 + 1
637 RETURN
638 END IF
639 CALL sscal( nru, one/nrmu,
640 $ z( irowu,icolz+i ), 2 )
641 IF( nrmu.NE.one .AND.
642 $ abs( nrmu-ortol )*sqrt2.GT.one )
643 $ THEN
644 DO j = 0, i-1
645 zjtji = -
sdot( nru, z( irowu, icolz+j ),
646 $ 2, z( irowu, icolz+i ), 2 )
647 CALL saxpy( nru, zjtji,
648 $ z( irowu, icolz+j ), 2,
649 $ z( irowu, icolz+i ), 2 )
650 END DO
651 nrmu =
snrm2( nru, z( irowu, icolz+i ), 2 )
652 CALL sscal( nru, one/nrmu,
653 $ z( irowu,icolz+i ), 2 )
654 END IF
655 END DO
656 DO i = 0, min( nsl-1, nrv-1 )
657 nrmv =
snrm2( nrv, z( irowv, icolz+i ), 2 )
658 IF( nrmv.EQ.zero ) THEN
659 info = n*2 + 1
660 RETURN
661 END IF
662 CALL sscal( nrv, -one/nrmv,
663 $ z( irowv,icolz+i ), 2 )
664 IF( nrmv.NE.one .AND.
665 $ abs( nrmv-ortol )*sqrt2.GT.one )
666 $ THEN
667 DO j = 0, i-1
668 zjtji = -
sdot( nrv, z( irowv, icolz+j ),
669 $ 2, z( irowv, icolz+i ), 2 )
670 CALL saxpy( nru, zjtji,
671 $ z( irowv, icolz+j ), 2,
672 $ z( irowv, icolz+i ), 2 )
673 END DO
674 nrmv =
snrm2( nrv, z( irowv, icolz+i ), 2 )
675 CALL sscal( nrv, one/nrmv,
676 $ z( irowv,icolz+i ), 2 )
677 END IF
678 END DO
679 IF( vutgk.EQ.zero .AND.
680 $ idptr.LT.idend .AND.
681 $ mod(ntgk,2).GT.0 ) THEN
682
683
684
685
686
687
688 split = .true.
689 z( irowz:irowz+ntgk-1,n+1 ) =
690 $ z( irowz:irowz+ntgk-1,ns+nsl )
691 z( irowz:irowz+ntgk-1,ns+nsl ) =
692 $ zero
693 END IF
694 END IF
695
696 nsl = min( nsl, nru )
697 sveq0 = .false.
698
699
700
701 DO i = 0, nsl-1
702 s( isbeg+i ) = abs( s( isbeg+i ) )
703 END DO
704
705
706
707 isbeg = isbeg + nsl
708 irowz = irowz + ntgk
709 icolz = icolz + nsl
710 irowu = irowz
711 irowv = irowz + 1
712 isplt = idptr + 1
713 ns = ns + nsl
714 nru = 0
715 nrv = 0
716 END IF
717 IF( irowz.LT.n*2 .AND. wantz ) THEN
718 z( 1:irowz-1, icolz ) = zero
719 END IF
720 END DO
721 IF( split .AND. wantz ) THEN
722
723
724
725
726 z( idbeg:idend-ntgk+1,isbeg-1 ) =
727 $ z( idbeg:idend-ntgk+1,isbeg-1 ) +
728 $ z( idbeg:idend-ntgk+1,n+1 )
729 z( idbeg:idend-ntgk+1,n+1 ) = 0
730 END IF
731 irowv = irowv - 1
732 irowu = irowu + 1
733 idbeg = ieptr + 1
734 sveq0 = .false.
735 split = .false.
736 END IF
737 END DO
738
739
740
741
742 DO i = 1, ns-1
743 k = 1
744 smin = s( 1 )
745 DO j = 2, ns + 1 - i
746 IF( s( j ).LE.smin ) THEN
747 k = j
748 smin = s( j )
749 END IF
750 END DO
751 IF( k.NE.ns+1-i ) THEN
752 s( k ) = s( ns+1-i )
753 s( ns+1-i ) = smin
754 IF( wantz )
CALL sswap( n*2, z( 1,k ), 1, z( 1,ns+1-i ), 1 )
755 END IF
756 END DO
757
758
759
760 IF( indsv ) THEN
761 k = iu - il + 1
762 IF( k.LT.ns ) THEN
763 s( k+1:ns ) = zero
764 IF( wantz ) z( 1:n*2,k+1:ns ) = zero
765 ns = k
766 END IF
767 END IF
768
769
770
771
772 IF( wantz ) THEN
773 DO i = 1, ns
774 CALL scopy( n*2, z( 1,i ), 1, work, 1 )
775 IF( lower ) THEN
776 CALL scopy( n, work( 2 ), 2, z( n+1,i ), 1 )
777 CALL scopy( n, work( 1 ), 2, z( 1 ,i ), 1 )
778 ELSE
779 CALL scopy( n, work( 2 ), 2, z( 1 ,i ), 1 )
780 CALL scopy( n, work( 1 ), 2, z( n+1,i ), 1 )
781 END IF
782 END DO
783 END IF
784
785 RETURN
786
787
788
subroutine xerbla(srname, info)
subroutine saxpy(n, sa, sx, incx, sy, incy)
SAXPY
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
real function sdot(n, sx, incx, sy, incy)
SDOT
integer function isamax(n, sx, incx)
ISAMAX
real function slamch(cmach)
SLAMCH
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.
logical function lsame(ca, cb)
LSAME
real(wp) function snrm2(n, x, incx)
SNRM2
subroutine sscal(n, sa, sx, incx)
SSCAL
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 matrice...
subroutine sswap(n, sx, incx, sy, incy)
SSWAP