3 SUBROUTINE pcgsepsubtst( WKNOWN, IBTYPE, JOBZ, RANGE, UPLO, N, VL,
4 $ VU, IL, IU, THRESH, ABSTOL, A, COPYA, B,
5 $ COPYB, Z, IA, JA, DESCA, WIN, WNEW,
6 $ IFAIL, ICLUSTR, GAP, IPREPAD, IPOSTPAD,
7 $ WORK, LWORK, RWORK, LRWORK, LWORK1,
8 $ IWORK, LIWORK, RESULT, TSTNRM, QTQNRM,
18 CHARACTER JOBZ, RANGE, UPLO
19 INTEGER IA, IBTYPE, IL, IPOSTPAD, IPREPAD, IU, JA,
20 $ LIWORK, LRWORK, LWORK, LWORK1, N, NOUT, RESULT
21 REAL ABSTOL, QTQNRM, THRESH, TSTNRM, VL, VU
24 INTEGER DESCA( * ), ICLUSTR( * ), IFAIL( * ),
26 REAL GAP( * ), RWORK( * ), WIN( * ), WNEW( * )
27 COMPLEX A( * ), B( * ), COPYA( * ), COPYB( * ),
223 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
224 $ MB_, NB_, RSRC_, CSRC_, LLD_
225 PARAMETER ( BLOCK_CYCLIC_2D = 1, dlen_ = 9, dtype_ = 1,
226 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
227 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
228 REAL PADVAL, FIVE, NEGONE
229 parameter( padval = 13.5285e+0, five = 5.0e+0,
232 PARAMETER ( CPADVAL = ( 13.989e+0, 1.93e+0 ) )
234 parameter( ipadval = 927 )
237 LOGICAL MISSLARGEST, MISSSMALLEST
238 INTEGER I, IAM, INDIWRK, INFO, ISIZEHEEVX, ISIZESUBTST,
239 $ isizetst, j, m, maxeigs, maxil, maxiu, maxsize,
240 $ minil, mq, mycol, myil, myrow, nclusters, np,
241 $ npcol, nprow, nq, nz, oldil, oldiu, oldnz, res,
242 $ rsizechk, rsizeheevx, rsizeqtq, rsizesubtst,
243 $ rsizetst, sizeheevx, sizemqrleft, sizemqrright,
244 $ sizeqrf, sizesubtst, sizetms, sizetst, valsize,
246 REAL EPS, ERROR, MAXERROR, MAXVU, MINERROR, MINVL,
247 $ NORMWIN, OLDVL, OLDVU, ORFAC, SAFMIN
250 INTEGER DESCZ( DLEN_ ), DSEED( 4 ), ITMP( 2 )
257 EXTERNAL LSAME, NUMROC, PSLAMCH
260 EXTERNAL blacs_gridinfo, clacpy,
descinit, igamn2d,
267 INTRINSIC abs,
max,
min, mod
271 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
273 CALL pclasizegsep( desca, iprepad, ipostpad, sizemqrleft,
274 $ sizemqrright, sizeqrf, sizetms, rsizeqtq,
275 $ rsizechk, sizeheevx, rsizeheevx, isizeheevx,
276 $ sizesubtst, rsizesubtst, isizesubtst, sizetst,
277 $ rsizetst, isizetst )
281 eps = pslamch( desca( ctxt_ ),
'Eps' )
282 safmin = pslamch( desca( ctxt_ ),
'Safe min' )
284 normwin = safmin / eps
286 $ normwin =
max( abs( win( 1 ) ), abs( win( n ) ), normwin )
297 DO 10 i = 1, lwork1, 1
298 rwork( i+iprepad ) = 14.3e+0
300 DO 20 i = 1, liwork, 1
301 iwork( i+iprepad ) = 14
303 DO 30 i = 1, lwork, 1
304 work( i+iprepad ) = ( 15.63e+0, 1.1e+0 )
308 wnew( i+iprepad ) = 3.14159e+0
311 iclustr( 1+iprepad ) = 139
313 IF( lsame( jobz,
'N' ) )
THEN
316 IF( lsame( range,
'A' ) )
THEN
318 ELSE IF( lsame( range,
'I' ) )
THEN
319 maxeigs = iu - il + 1
321 minvl = vl - normwin*five*eps - abstol
322 maxvu = vu + normwin*five*eps + abstol
326 IF( win( i ).LT.minvl )
328 IF( win( i ).LE.maxvu )
332 maxeigs = maxiu - minil + 1
337 CALL descinit( descz, desca( m_ ), desca( n_ ), desca( mb_ ),
338 $ desca( nb_ ), desca( rsrc_ ), desca( csrc_ ),
339 $ desca( ctxt_ ), desca( lld_ ), info )
341 CALL blacs_gridinfo( desca( ctxt_ ), nprow, npcol, myrow, mycol )
342 indiwrk = 1 + iprepad + nprow*npcol + 1
345 IF( myrow.EQ.0 .AND. mycol.EQ.0 )
351 IF( myrow.GE.nprow .OR. myrow.LT.0 )
362 $ dseed, win, maxsize, vecsize, valsize )
364 np = numroc( n, desca( mb_ ), myrow, 0, nprow )
365 nq = numroc( n, desca( nb_ ), mycol, 0, npcol )
366 mq = numroc( maxeigs, desca( nb_ ), mycol, 0, npcol )
368 CALL clacpy(
'A', np, nq, copya, desca( lld_ ), a( 1+iprepad ),
371 CALL clacpy(
'A', np, nq, copyb, desca( lld_ ), b( 1+iprepad ),
374 CALL pcfillpad( desca( ctxt_ ), np, nq, b, desca( lld_ ), iprepad,
375 $ ipostpad, cpadval+1.0e+2 )
377 CALL pcfillpad( desca( ctxt_ ), np, nq, a, desca( lld_ ), iprepad,
378 $ ipostpad, cpadval )
380 CALL pcfillpad( descz( ctxt_ ), np, mq, z, descz( lld_ ), iprepad,
381 $ ipostpad, cpadval+1.0e+0 )
383 CALL psfillpad( desca( ctxt_ ), n, 1, wnew, n, iprepad, ipostpad,
386 CALL psfillpad( desca( ctxt_ ), nprow*npcol, 1, gap, nprow*npcol,
387 $ iprepad, ipostpad, padval+3.0e+0 )
389 CALL psfillpad( desca( ctxt_ ), lwork1, 1, rwork, lwork1, iprepad,
390 $ ipostpad, padval+4.0e+0 )
392 CALL pifillpad( desca( ctxt_ ), liwork, 1, iwork, liwork, iprepad,
393 $ ipostpad, ipadval )
395 CALL pifillpad( desca( ctxt_ ), n, 1, ifail, n, iprepad, ipostpad,
398 CALL pifillpad( desca( ctxt_ ), 2*nprow*npcol, 1, iclustr,
399 $ 2*nprow*npcol, iprepad, ipostpad, ipadval )
401 CALL pcfillpad( desca( ctxt_ ), lwork, 1, work, lwork, iprepad,
402 $ ipostpad, cpadval+4.1e+0 )
408 DO 60 j = 1, maxeigs, 1
409 CALL pcelset( z( 1+iprepad ), i, j, desca,
410 $ ( 13.0e+0, 1.34e+0 ) )
419 CALL pchegvx( ibtype, jobz, range, uplo, n, a( 1+iprepad ), ia,
420 $ ja, desca, b( 1+iprepad ), ia, ja, desca, vl, vu,
421 $ il, iu, abstol, m, nz, wnew( 1+iprepad ), orfac,
422 $ z( 1+iprepad ), ia, ja, desca, work( 1+iprepad ),
423 $ sizeheevx, rwork( 1+iprepad ), lwork1,
424 $ iwork( 1+iprepad ), liwork, ifail( 1+iprepad ),
425 $ iclustr( 1+iprepad ), gap( 1+iprepad ), info )
429 IF( thresh.LE.0 )
THEN
432 CALL pcchekpad( desca( ctxt_ ),
'PCHEGVX-B', np, nq, b,
433 $ desca( lld_ ), iprepad, ipostpad,
436 CALL pcchekpad( desca( ctxt_ ),
'PCHEGVX-A', np, nq, a,
437 $ desca( lld_ ), iprepad, ipostpad, cpadval )
439 CALL pcchekpad( descz( ctxt_ ),
'PCHEGVX-Z', np, mq, z,
440 $ descz( lld_ ), iprepad, ipostpad,
443 CALL pschekpad( desca( ctxt_ ),
'PCHEGVX-WNEW', n, 1, wnew, n,
444 $ iprepad, ipostpad, padval+2.0e+0 )
446 CALL pschekpad( desca( ctxt_ ),
'PCHEGVX-GAP', nprow*npcol, 1,
447 $ gap, nprow*npcol, iprepad, ipostpad,
450 CALL pschekpad( desca( ctxt_ ),
'PCHEGVX-rWORK', lwork1, 1,
451 $ rwork, lwork1, iprepad, ipostpad,
454 CALL pcchekpad( desca( ctxt_ ),
'PCHEGVX-WORK', lwork, 1, work,
455 $ lwork, iprepad, ipostpad, cpadval+4.1e+0 )
457 CALL pichekpad( desca( ctxt_ ),
'PCHEGVX-IWORK', liwork, 1,
458 $ iwork, liwork, iprepad, ipostpad, ipadval )
460 CALL pichekpad( desca( ctxt_ ),
'PCHEGVX-IFAIL', n, 1, ifail,
461 $ n, iprepad, ipostpad, ipadval )
463 CALL pichekpad( desca( ctxt_ ),
'PCHEGVX-ICLUSTR',
464 $ 2*nprow*npcol, 1, iclustr, 2*nprow*npcol,
465 $ iprepad, ipostpad, ipadval )
470 IF( lsame( range,
'A' ) )
THEN
472 $ dseed, wnew( 1+iprepad ), maxsize,
485 CALL igamn2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp, 1, 1, 1,
487 CALL igamx2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp( 2 ), 1, 1,
491 IF( itmp( 1 ).NE.itmp( 2 ) )
THEN
493 $
WRITE( nout, fmt = * )
494 $
'Different processes return different INFO'
496 ELSE IF( mod( info, 2 ).EQ.1 .OR. info.GT.7 .OR. info.LT.0 )
499 $
WRITE( nout, fmt = 9999 )info
501 ELSE IF( mod( info / 2, 2 ).EQ.1 .AND. lwork1.GE.maxsize )
THEN
503 $
WRITE( nout, fmt = 9996 )info
505 ELSE IF( mod( info / 4, 2 ).EQ.1 .AND. lwork1.GE.vecsize )
THEN
507 $
WRITE( nout, fmt = 9996 )info
512 IF( lsame( jobz,
'V' ) .AND. ( iclustr( 1+iprepad ).NE.
513 $ 0 ) .AND. ( mod( info / 2, 2 ).NE.1 ) )
THEN
515 $
WRITE( nout, fmt = 9995 )
521 IF( ( m.LT.0 ) .OR. ( m.GT.n ) )
THEN
523 $
WRITE( nout, fmt = 9994 )
525 ELSE IF( lsame( range,
'A' ) .AND. ( m.NE.n ) )
THEN
527 $
WRITE( nout, fmt = 9993 )
529 ELSE IF( lsame( range,
'I' ) .AND. ( m.NE.iu-il+1 ) )
THEN
531 $
WRITE( nout, fmt = 9992 )
533 ELSE IF( lsame( jobz,
'V' ) .AND.
534 $ ( .NOT.( lsame( range,
'V' ) ) ) .AND. ( m.NE.nz ) )
537 $
WRITE( nout, fmt = 9991 )
543 IF( lsame( jobz,
'V' ) )
THEN
544 IF( lsame( range,
'V' ) )
THEN
547 $
WRITE( nout, fmt = 9990 )
550 IF( nz.LT.m .AND. mod( info / 4, 2 ).NE.1 )
THEN
552 $
WRITE( nout, fmt = 9989 )
558 $
WRITE( nout, fmt = 9988 )
563 IF( result.EQ.0 )
THEN
570 CALL igamn2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp, 1, 1, 1,
572 CALL igamx2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp( 2 ), 1,
575 IF( itmp( 1 ).NE.itmp( 2 ) )
THEN
577 $
WRITE( nout, fmt = 9987 )
584 rwork( i ) = wnew( i+iprepad )
585 rwork( i+m ) = wnew( i+iprepad )
588 CALL sgamn2d( desca( ctxt_ ),
'a',
' ', m, 1, rwork, m,
590 CALL sgamx2d( desca( ctxt_ ),
'a',
' ', m, 1,
591 $ rwork( 1+m ), m, 1, 1, -1, -1, 0 )
595 IF( result.EQ.0 .AND. ( abs( rwork( i )-rwork( m+
596 $ i ) ).GT.five*eps*abs( rwork( i ) ) ) )
THEN
598 $
WRITE( nout, fmt = 9986 )
607 IF( lsame( jobz,
'V' ) )
THEN
609 DO 100 i = 0, nprow*npcol - 1
610 IF( iclustr( 1+iprepad+2*i ).EQ.0 )
612 nclusters = nclusters + 1
615 itmp( 1 ) = nclusters
616 itmp( 2 ) = nclusters
618 CALL igamn2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp, 1, 1, 1,
620 CALL igamx2d( desca( ctxt_ ),
'a',
' ', 1, 1, itmp( 2 ), 1,
623 IF( itmp( 1 ).NE.itmp( 2 ) )
THEN
625 $
WRITE( nout, fmt = 9985 )
631 DO 120 i = 1, nclusters
632 iwork( indiwrk+i ) = iclustr( i+iprepad )
633 iwork( indiwrk+i+nclusters ) = iclustr( i+iprepad )
635 CALL igamn2d( desca( ctxt_ ),
'a',
' ', nclusters*2+1, 1,
636 $ iwork( indiwrk+1 ), nclusters*2+1, 1, 1,
638 CALL igamx2d( desca( ctxt_ ),
'a',
' ', nclusters*2+1, 1,
639 $ iwork( indiwrk+1+nclusters ),
640 $ nclusters*2+1, 1, 1, -1, -1, 0 )
643 DO 130 i = 1, nclusters
644 IF( result.EQ.0 .AND. iwork( indiwrk+i ).NE.
645 $ iwork( indiwrk+nclusters+i ) )
THEN
647 $
WRITE( nout, fmt = 9984 )
652 IF( iclustr( 1+iprepad+nclusters*2 ).NE.0 )
THEN
654 $
WRITE( nout, fmt = 9983 )
661 CALL igamx2d( desca( ctxt_ ),
'a',
' ', 1, 1, result, 1, 1, 1,
677 IF( lsame( jobz,
'V' ) )
THEN
681 CALL psfillpad( desca( ctxt_ ), rsizechk, 1, rwork,
682 $ rsizechk, iprepad, ipostpad, 4.3e+0 )
684 CALL pcgsepchk( ibtype, n, nz, copya, ia, ja, desca, copyb,
685 $ ia, ja, desca, thresh, z( 1+iprepad ), ia,
686 $ ja, descz, a( 1+iprepad ), ia, ja, desca,
687 $ wnew( 1+iprepad ), rwork( 1+iprepad ),
688 $ rsizechk, tstnrm, res )
690 CALL pschekpad( desca( ctxt_ ),
'PCGSEPCHK-rWORK', rsizechk,
691 $ 1, rwork, rsizechk, iprepad, ipostpad,
706 IF( lsame( range,
'V' ) )
THEN
711 IF( lsame( range,
'A' ) )
THEN
723 DO 150 myil = minil, maxil
728 misssmallest = .true.
729 IF( .NOT.lsame( range,
'V' ) .OR. ( myil.EQ.1 ) )
730 $ misssmallest = .false.
731 IF( misssmallest .AND. ( win( myil-1 ).LT.vl+normwin*
732 $ five*thresh*eps ) )misssmallest = .false.
734 IF( .NOT.lsame( range,
'V' ) .OR. ( myil.EQ.maxil ) )
735 $ misslargest = .false.
736 IF( misslargest .AND. ( win( myil+m ).GT.vu-normwin*five*
737 $ thresh*eps ) )misslargest = .false.
738 IF( .NOT.misssmallest )
THEN
739 IF( .NOT.misslargest )
THEN
744 error = abs( win( i+myil-1 )-wnew( i+iprepad ) )
745 maxerror =
max( maxerror, error )
748 minerror =
min( maxerror, minerror )
759 IF( lsame( jobz,
'V' ) .AND. lsame( range,
'A' ) )
THEN
760 IF( minerror.GT.normwin*five*five*thresh*eps )
THEN
762 $
WRITE( nout, fmt = 9997 )minerror, normwin
766 IF( minerror.GT.normwin*five*thresh*eps )
THEN
768 $
WRITE( nout, fmt = 9997 )minerror, normwin
777 IF( il.NE.oldil .OR. iu.NE.oldiu .OR. vl.NE.oldvl .OR. vu.NE.
780 $
WRITE( nout, fmt = 9982 )
784 IF( lsame( jobz,
'N' ) .AND. ( nz.NE.oldnz ) )
THEN
786 $
WRITE( nout, fmt = 9981 )
794 CALL igamx2d( desca( ctxt_ ),
'a',
' ', 1, 1, result, 1, 1, 1, -1,
802 9999
FORMAT(
'PCHEGVX returned INFO=', i7 )
803 9998
FORMAT(
'PCSEPQTQ returned INFO=', i7 )
804 9997
FORMAT(
'PCGSEPSUBTST minerror =', d11.2,
' normwin=', d11.2 )
805 9996
FORMAT(
'PCHEGVX returned INFO=', i7,
806 $
' despite adequate workspace' )
807 9995
FORMAT( .NE..NE.
'ICLUSTR(1)0 but mod(INFO/2,2)1' )
808 9994
FORMAT(
'M not in the range 0 to N' )
809 9993
FORMAT(
'M not equal to N' )
810 9992
FORMAT(
'M not equal to IU-IL+1' )
811 9991
FORMAT(
'M not equal to NZ' )
812 9990
FORMAT(
'NZ > M' )
813 9989
FORMAT(
'NZ < M' )
814 9988
FORMAT(
'NZ not equal to M' )
815 9987
FORMAT(
'Different processes return different values for M' )
816 9986
FORMAT(
'Different processes return different eigenvalues' )
817 9985
FORMAT(
'Different processes return ',
818 $
'different numbers of clusters' )
819 9984
FORMAT(
'Different processes return different clusters' )
820 9983
FORMAT(
'ICLUSTR not zero terminated' )
821 9982
FORMAT(
'IL, IU, VL or VU altered by PCHEGVX' )
822 9981
FORMAT(
'NZ altered by PCHEGVX with JOBZ=N' )