1 SUBROUTINE pdsyevr( JOBZ, RANGE, UPLO, N, A, IA, JA,
2 $ DESCA, VL, VU, IL, IU, M, NZ, W, Z, IZ,
3 $ JZ, DESCZ, WORK, LWORK, IWORK, LIWORK,
13 CHARACTER JOBZ, RANGE, UPLO
14 INTEGER IA, IL, INFO, IU, IZ, JA, JZ, LIWORK, LWORK, M,
16 DOUBLE PRECISION VL, VU
19 INTEGER DESCA( * ), DESCZ( * ), IWORK( * )
20 DOUBLE PRECISION A( * ), W( * ), WORK( * ), Z( * )
294 INTEGER CTXT_, M_, N_,
295 $ MB_, NB_, RSRC_, CSRC_
296 PARAMETER ( CTXT_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
297 $ rsrc_ = 7, csrc_ = 8 )
298 DOUBLE PRECISION ZERO
299 parameter( zero = 0.0d0 )
302 LOGICAL ALLEIG, COLBRT, DOBCST, FINISH, FIRST, INDEIG,
303 $ LOWER, LQUERY, VALEIG, VSTART, WANTZ
304 INTEGER ANB, DOL, DOU, DSTCOL, DSTROW, EIGCNT, FRSTCL,
305 $ i, iarow, ictxt, iil, iinderr, iindwlc, iinfo,
306 $ iiu, im, indd, indd2, inde, inde2, inderr,
307 $ indilu, indrw, indtau, indwlc, indwork, ipil,
308 $ ipiu, iproc, izrow, lastcl, lengthi, lengthi2,
309 $ liwmin, llwork, lwmin, lwopt, maxcls, mq00,
310 $ mycol, myil, myiu, myproc, myrow, mz, nb,
311 $ ndepth, needil, neediu, nnp, np00, npcol,
312 $ nprocs, nprow, nps, nsplit, nsytrd_lwopt,
313 $ offset, parity, rlengthi, rlengthi2, rstarti,
314 $ size1, size2, sqnpc, srccol, srcrow, starti,
317 DOUBLE PRECISION PIVMIN, SAFMIN, SCALE, VLL, VUU, WL,
321 INTEGER IDUM1( 4 ), IDUM2( 4 )
325 INTEGER ICEIL, INDXG2P, NUMROC, PJLAENV
326 DOUBLE PRECISION PDLAMCH
327 EXTERNAL iceil, indxg2p, lsame, numroc, pdlamch,
331 EXTERNAL blacs_gridinfo,
chk1mat, dcopy, dgebr2d,
332 $ dgebs2d, dgerv2d, dgesd2d, dlarrc,
dlasrt2,
334 $ igebs2d, igerv2d, igesd2d, igsum2d,
pchk1mat,
339 INTRINSIC abs, dble, ichar, int,
max,
min, mod, sqrt
351 wantz = lsame( jobz,
'V' )
352 lower = lsame( uplo,
'L' )
353 alleig = lsame( range,
'A' )
354 valeig = lsame( range,
'V' )
355 indeig = lsame( range,
'I' )
356 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
363 ictxt = desca( ctxt_ )
364 safmin = pdlamch( ictxt,
'Safe minimum' )
377 llwork = lwork - indwork + 1
384 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
387 nprocs = nprow * npcol
388 myproc = myrow * npcol + mycol
389 IF( nprow.EQ.-1 )
THEN
390 info = -( 800+ctxt_ )
391 ELSE IF( wantz )
THEN
392 IF( ictxt.NE.descz( ctxt_ ) )
THEN
393 info = -( 2100+ctxt_ )
404 ELSE IF ( indeig )
THEN
413 np00 = numroc( n, nb, 0, 0, nprow )
414 mq00 = numroc( mz, nb, 0, 0, npcol )
415 indrw = indwork +
max(18*n, np00*mq00 + 2*nb*nb)
416 lwmin = indrw - 1 + (iceil(mz, nprocs) + 2)*n
418 indrw = indwork + 12*n
422 lwmin =
max(3, lwmin)
424 anb = pjlaenv( ictxt, 3,
'PDSYTTRD',
'L', 0, 0, 0, 0 )
425 sqnpc = int( sqrt( dble( nprocs ) ) )
426 nps =
max( numroc( n, 1, 0, 0, sqnpc ), 2*anb )
427 nsytrd_lwopt = 2*( anb+1 )*( 4*nps+2 ) + ( nps+4 )*nps
428 lwopt =
max( lwopt, 5*n+nsytrd_lwopt )
430 size1 = indrw - indwork
437 nnp =
max( n, nprocs+1, 4 )
439 liwmin = 12*nnp + 2*n
441 liwmin = 10*nnp + 2*n
450 indilu = liwmin - 2*nprocs + 1
460 CALL chk1mat( n, 4, n, 4, ia, ja, desca, 8, info )
462 $
CALL chk1mat( n, 4, n, 4, iz, jz, descz, 21, info )
465 IF( .NOT.( wantz .OR. lsame( jobz,
'N' ) ) )
THEN
467 ELSE IF( .NOT.( alleig .OR. valeig .OR. indeig ) )
THEN
469 ELSE IF( .NOT.( lower .OR. lsame( uplo,
'U' ) ) )
THEN
471 ELSE IF( mod( ia-1, desca( mb_ ) ).NE.0 )
THEN
473 ELSE IF( valeig .AND. n.GT.0 .AND. vu.LE.vl )
THEN
475 ELSE IF( indeig .AND. ( il.LT.1 .OR. il.GT.
max( 1, n ) ) )
478 ELSE IF( indeig .AND. ( iu.LT.
min( n, il ) .OR. iu.GT.n ) )
481 ELSE IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
483 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery )
THEN
485 ELSE IF( desca( mb_ ).NE.desca( nb_ ) )
THEN
489 iarow = indxg2p( 1, desca( nb_ ), myrow,
490 $ desca( rsrc_ ), nprow )
491 izrow = indxg2p( 1, desca( nb_ ), myrow,
492 $ descz( rsrc_ ), nprow )
493 IF( iarow.NE.izrow )
THEN
495 ELSE IF( mod( ia-1, desca( mb_ ) ).NE.
496 $ mod( iz-1, descz( mb_ ) ) )
THEN
498 ELSE IF( desca( m_ ).NE.descz( m_ ) )
THEN
500 ELSE IF( desca( n_ ).NE.descz( n_ ) )
THEN
502 ELSE IF( desca( mb_ ).NE.descz( mb_ ) )
THEN
504 ELSE IF( desca( nb_ ).NE.descz( nb_ ) )
THEN
506 ELSE IF( desca( rsrc_ ).NE.descz( rsrc_ ) )
THEN
507 info = -( 2100+rsrc_ )
508 ELSE IF( desca( csrc_ ).NE.descz( csrc_ ) )
THEN
509 info = -( 2100+csrc_ )
510 ELSE IF( ictxt.NE.descz( ctxt_ ) )
THEN
511 info = -( 2100+ctxt_ )
517 idum1( 2 ) = ichar(
'L' )
519 idum1( 2 ) = ichar(
'U' )
523 idum1( 3 ) = ichar(
'A' )
524 ELSE IF( indeig )
THEN
525 idum1( 3 ) = ichar(
'I' )
527 idum1( 3 ) = ichar(
'V' )
537 idum1( 1 ) = ichar(
'V' )
538 CALL pchk2mat( n, 4, n, 4, ia, ja, desca, 8, n, 4, n, 4, iz,
539 $ jz, descz, 21, 4, idum1, idum2, info )
541 idum1( 1 ) = ichar(
'N' )
542 CALL pchk1mat( n, 4, n, 4, ia, ja, desca, 8, 4, idum1,
545 work( 1 ) = dble( lwopt )
550 CALL pxerbla( ictxt,
'PDSYEVR', -info )
552 ELSE IF( lquery )
THEN
566 work( 1 ) = dble( lwopt )
589 CALL pdsyntrd( uplo, n, a, ia, ja, desca, work( indd ),
590 $ work( inde ), work( indtau ), work( indwork ),
594 IF (iinfo .NE. 0)
THEN
595 CALL pxerbla( ictxt,
'PDSYNTRD', -iinfo )
605 IF( ia.EQ.1 .AND. ja.EQ.1 .AND.
606 $ desca( rsrc_ ).EQ.0 .AND. desca( csrc_ ).EQ.0 )
608 CALL pdlared1d( n, ia, ja, desca, work( indd ), work( indd2 ),
609 $ work( indwork ), llwork )
611 CALL pdlared1d( n, ia, ja, desca, work( inde ), work( inde2 ),
612 $ work( indwork ), llwork )
617 CALL pdelget(
'A',
' ', work( indd2+i-1 ), a, i+ia-1,
620 IF( lsame( uplo,
'U' ) )
THEN
622 CALL pdelget(
'A',
' ', work( inde2+i-1 ), a, i+ia-1,
627 CALL pdelget(
'A',
' ', work( inde2+i-1 ), a, i+ia,
644 ELSE IF ( indeig )
THEN
647 ELSE IF ( valeig )
THEN
648 CALL dlarrc(
'T', n, vll, vuu, work( indd2 ),
649 $ work( inde2 + offset ), safmin, eigcnt, iil, iiu, info)
660 work( 1 ) = dble( lwopt )
678 CALL pmpim2( iil, iiu, nprocs,
679 $ iwork(indilu), iwork(indilu+nprocs) )
683 myil = iwork(indilu+myproc)
684 myiu = iwork(indilu+nprocs+myproc)
687 zoffset =
max(0, myil - iil - 1)
688 first = ( myil .EQ. iil )
701 IF ( myil.GT.0 )
THEN
703 dou = myiu - myil + 1
704 CALL dstegr2( jobz,
'I', n, work( indd2 ),
705 $ work( inde2+offset ), vll, vuu, myil, myiu,
706 $ im, w( 1 ), work( indrw ), n,
708 $ iwork( 1 ), work( indwork ), size1,
709 $ iwork( 2*n+1 ), size2,
710 $ dol, dou, zoffset, iinfo )
715 w( myil-iil+i ) = w( i )
720 IF (iinfo .NE. 0)
THEN
721 CALL pxerbla( ictxt,
'DSTEGR2', -iinfo )
724 ELSEIF ( wantz .AND. nprocs.EQ.1 )
THEN
729 IF ( myil.GT.0 )
THEN
732 CALL dstegr2( jobz,
'I', n, work( indd2 ),
733 $ work( inde2+offset ), vll, vuu, iil, iiu,
734 $ im, w( 1 ), work( indrw ), n,
736 $ iwork( 1 ), work( indwork ), size1,
737 $ iwork( 2*n+1 ), size2, dol, dou,
740 IF (iinfo .NE. 0)
THEN
741 CALL pxerbla( ictxt,
'DSTEGR2', -iinfo )
744 ELSEIF ( wantz )
THEN
752 IF ( myil.GT.0 )
THEN
755 CALL dstegr2a( jobz,
'I', n, work( indd2 ),
756 $ work( inde2+offset ), vll, vuu, iil, iiu,
757 $ im, w( 1 ), work( indrw ), n,
758 $ n, work( indwork ), size1,
759 $ iwork( 2*n+1 ), size2, dol,
760 $ dou, needil, neediu,
761 $ inderr, nsplit, pivmin, scale, wl, wu,
764 IF (iinfo .NE. 0)
THEN
765 CALL pxerbla( ictxt,
'DSTEGR2A', -iinfo )
777 iinderr = indwork + inderr - 1
795 IF (myproc .EQ. (i - 1))
THEN
801 lengthi = myiu - myil + 1
806 CALL igesd2d( ictxt, 2, 1, iwork, 2,
808 IF (( starti.GE.1 ) .AND. ( lengthi.GE.1 ))
THEN
811 CALL dcopy(lengthi,w( starti ),1,
814 CALL dcopy(lengthi,work( iinderr+starti-1 ),1,
815 $ work( indd+lengthi ), 1)
817 CALL dgesd2d( ictxt, lengthi2,
818 $ 1, work( indd ), lengthi2,
821 ELSE IF (myproc .EQ. 0)
THEN
822 srcrow = (i-1) / npcol
823 srccol = mod(i-1, npcol)
824 CALL igerv2d( ictxt, 2, 1, iwork, 2,
828 IF (( starti.GE.1 ) .AND. ( lengthi.GE.1 ))
THEN
831 CALL dgerv2d( ictxt, lengthi2, 1,
832 $ work(indd), lengthi2, srcrow, srccol )
834 CALL dcopy( lengthi, work(indd), 1,
837 CALL dcopy(lengthi,work(indd+lengthi),1,
838 $ work( iinderr+starti-1 ), 1)
842 lengthi = iiu - iil + 1
843 lengthi2 = lengthi * 2
844 IF (myproc .EQ. 0)
THEN
846 CALL dcopy(lengthi,w ,1, work( indd ), 1)
847 CALL dcopy(lengthi,work( iinderr ),1,
848 $ work( indd+lengthi ), 1)
849 CALL dgebs2d( ictxt,
'A',
' ', lengthi2, 1,
850 $ work(indd), lengthi2 )
854 CALL dgebr2d( ictxt,
'A',
' ', lengthi2, 1,
855 $ work(indd), lengthi2, srcrow, srccol )
856 CALL dcopy( lengthi, work(indd), 1, w, 1)
857 CALL dcopy(lengthi,work(indd+lengthi),1,
858 $ work( iinderr ), 1)
865 IF( (nprocs.GT.1).AND.(myil.GT.0) )
THEN
866 CALL pmpcol( myproc, nprocs, iil, needil, neediu,
867 $ iwork(indilu), iwork(indilu+nprocs),
868 $ colbrt, frstcl, lastcl )
876 DO 47 iproc = frstcl, lastcl
877 IF (myproc .EQ. iproc)
THEN
880 lengthi = myiu - myil + 1
883 IF ((starti.GE.1) .AND. (lengthi.GE.1))
THEN
885 CALL dcopy(lengthi,w( starti ),1,
889 $ work( iinderr+starti-1 ),1,
890 $ work(indd+lengthi), 1)
893 DO 46 i = frstcl, lastcl
894 IF(i.EQ.myproc)
GOTO 46
896 dstcol = mod(i, npcol)
897 CALL igesd2d( ictxt, 2, 1, iwork, 2,
899 IF ((starti.GE.1) .AND. (lengthi.GE.1))
THEN
902 CALL dgesd2d( ictxt, lengthi2,
903 $ 1, work(indd), lengthi2,
908 srcrow = iproc / npcol
909 srccol = mod(iproc, npcol)
910 CALL igerv2d( ictxt, 2, 1, iwork, 2,
914 IF ((rstarti.GE.1 ) .AND. (rlengthi.GE.1 ))
THEN
915 rlengthi2 = 2*rlengthi
916 CALL dgerv2d( ictxt, rlengthi2, 1,
917 $ work(inde), rlengthi2,
920 CALL dcopy( rlengthi, work(inde), 1,
923 CALL dcopy(rlengthi,work(inde+rlengthi),1,
924 $ work( iinderr+rstarti-1 ), 1)
939 IF ( myil.GT.0 )
THEN
940 CALL dstegr2b( jobz, n, work( indd2 ),
941 $ work( inde2+offset ),
942 $ im, w( 1 ), work( indrw ), n, n,
943 $ iwork( 1 ), work( indwork ), size1,
944 $ iwork( 2*n+1 ), size2, dol,
945 $ dou, needil, neediu, indwlc,
946 $ pivmin, scale, wl, wu,
948 $ maxcls, ndepth, parity, zoffset, iinfo )
949 iindwlc = indwork + indwlc - 1
951 IF((needil.LT.dol).OR.(neediu.GT.dou))
THEN
952 CALL pmpcol( myproc, nprocs, iil, needil, neediu,
953 $ iwork(indilu), iwork(indilu+nprocs),
954 $ colbrt, frstcl, lastcl )
965 DO 147 iproc = frstcl, lastcl
966 IF (myproc .EQ. iproc)
THEN
970 lengthi = myiu - myil + 1
975 IF ((starti.GE.1).AND.(lengthi.GE.1))
THEN
978 $ work( iindwlc+starti-1 ),1,
982 $ work( iinderr+starti-1 ),1,
983 $ work(indd+lengthi), 1)
986 DO 146 i = frstcl, lastcl
987 IF(i.EQ.myproc)
GOTO 146
989 dstcol = mod(i, npcol)
990 CALL igesd2d( ictxt, 2, 1, iwork, 2,
992 IF ((starti.GE.1).AND.(lengthi.GE.1))
THEN
995 CALL dgesd2d( ictxt, lengthi2,
996 $ 1, work(indd), lengthi2,
1001 srcrow = iproc / npcol
1002 srccol = mod(iproc, npcol)
1003 CALL igerv2d( ictxt, 2, 1, iwork, 2,
1007 IF ((rstarti.GE.1).AND.(rlengthi.GE.1))
THEN
1008 rlengthi2 = 2*rlengthi
1009 CALL dgerv2d( ictxt,rlengthi2, 1,
1010 $ work(inde),rlengthi2,
1013 CALL dcopy(rlengthi, work(inde), 1,
1014 $ work( iindwlc+rstarti-1 ), 1)
1016 CALL dcopy(rlengthi,work(inde+rlengthi),1,
1017 $ work( iinderr+rstarti-1 ), 1)
1025 IF (iinfo .NE. 0)
THEN
1026 CALL pxerbla( ictxt,
'DSTEGR2B', -iinfo )
1047 IF (myproc .EQ. (i - 1))
THEN
1050 starti = myil - iil + 1
1053 lengthi = myiu - myil + 1
1058 CALL igesd2d( ictxt, 2, 1, iwork, 2,
1060 IF ((starti.GE.1).AND.(lengthi.GE.1))
THEN
1061 CALL dgesd2d( ictxt, lengthi,
1062 $ 1, w( starti ), lengthi,
1065 ELSE IF (myproc .EQ. 0)
THEN
1066 srcrow = (i-1) / npcol
1067 srccol = mod(i-1, npcol)
1068 CALL igerv2d( ictxt, 2, 1, iwork, 2,
1072 IF ((starti.GE.1).AND.(lengthi.GE.1))
THEN
1073 CALL dgerv2d( ictxt, lengthi, 1,
1074 $ w( starti ), lengthi, srcrow, srccol )
1081 CALL igsum2d( ictxt,
'A',
' ', 1, 1, m, 1, -1, -1 )
1084 IF (myproc .EQ. 0)
THEN
1086 CALL dgebs2d( ictxt,
'A',
' ', m, 1, w, m )
1090 CALL dgebr2d( ictxt,
'A',
' ', m, 1,
1091 $ w, m, srcrow, srccol )
1098 iwork( nprocs+1+i ) = i
1100 CALL dlasrt2(
'I', m, w, iwork( nprocs+2 ), iinfo )
1101 IF (iinfo.NE.0)
THEN
1102 CALL pxerbla( ictxt,
'DLASRT2', -iinfo )
1113 iwork( m+nprocs+1+iwork( nprocs+1+i ) ) = i
1117 DO 180 i = 1, nprocs
1120 ipil = iwork(indilu+i-1)
1121 ipiu = iwork(indilu+nprocs+i-1)
1122 IF (ipil .EQ. 0)
THEN
1123 iwork( i + 1 ) = iwork( i )
1125 iwork( i + 1 ) = iwork( i ) + ipiu - ipil + 1
1130 CALL pdlaevswp(n, work( indrw ), n, z, iz, jz,
1131 $ descz, iwork( 1 ), iwork( nprocs+m+2 ), work( indwork ),
1134 CALL pdlaevswp(n, work( indrw + n ), n, z, iz, jz,
1135 $ descz, iwork( 1 ), iwork( nprocs+m+2 ), work( indwork ),
1148 CALL pdormtr(
'L', uplo,
'N', n, nz, a, ia, ja, desca,
1149 $ work( indtau ), z, iz, jz, descz,
1150 $ work( indwork ), size1, iinfo )
1152 IF (iinfo.NE.0)
THEN
1153 CALL pxerbla( ictxt,
'PDORMTR', -iinfo )
1160 work( 1 ) = dble( lwopt )