299 SUBROUTINE sdrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
300 $ alphar, alphai, beta, vl, vr, ilo, ihi, lscale,
301 $ rscale, s, stru, dif, diftru, work, lwork,
302 $ iwork, liwork, result, bwork, info )
310 INTEGER ihi, ilo, info, lda, liwork, lwork, nin, nout,
317 REAL a( lda, * ), ai( lda, * ), alphai( * ),
318 $ alphar( * ), b( lda, * ), beta( * ),
319 $ bi( lda, * ), dif( * ), diftru( * ),
320 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
321 $ stru( * ), vl( lda, * ), vr( lda, * ),
328 REAL zero, one, ten, tnth, half
329 parameter( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
330 $ tnth = 1.0e-1, half = 0.5d+0 )
333 INTEGER i, iptype, iwa, iwb, iwx, iwy, j, linfo,
334 $ maxwrk, minwrk, n, nerrs, nmax, nptknt, ntestt
335 REAL abnorm, anorm, bnorm, ratio1, ratio2, thrsh2,
350 INTRINSIC abs, max, sqrt
360 IF( nsize.LT.0 )
THEN
362 ELSE IF( thresh.LT.zero )
THEN
364 ELSE IF( nin.LE.0 )
THEN
366 ELSE IF( nout.LE.0 )
THEN
368 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
370 ELSE IF( liwork.LT.nmax+6 )
THEN
382 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
383 minwrk = 2*nmax*nmax + 12*nmax + 16
384 maxwrk = 6*nmax + nmax*
ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
386 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
390 IF( lwork.LT.minwrk )
394 CALL
xerbla(
'SDRGVX', -info )
414 weight( 4 ) = one / weight( 2 )
415 weight( 5 ) = one / weight( 1 )
425 CALL
slatm6( iptype, 5, a, lda, b, vr, lda, vl,
426 $ lda, weight( iwa ), weight( iwb ),
427 $ weight( iwx ), weight( iwy ), stru,
434 CALL
slacpy(
'F', n, n, a, lda, ai, lda )
435 CALL
slacpy(
'F', n, n, b, lda, bi, lda )
437 CALL
sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
438 $ lda, alphar, alphai, beta, vl, lda,
439 $ vr, lda, ilo, ihi, lscale, rscale,
440 $ anorm, bnorm, s, dif, work, lwork,
441 $ iwork, bwork, linfo )
442 IF( linfo.NE.0 )
THEN
444 WRITE( nout, fmt = 9999 )
'SGGEVX', linfo, n,
451 CALL
slacpy(
'Full', n, n, ai, lda, work, n )
452 CALL
slacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
454 abnorm =
slange(
'Fro', n, 2*n, work, n, work )
459 CALL
sget52( .true., n, a, lda, b, lda, vl, lda,
460 $ alphar, alphai, beta, work,
462 IF( result( 2 ).GT.thresh )
THEN
463 WRITE( nout, fmt = 9998 )
'Left',
'SGGEVX',
464 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
468 CALL
sget52( .false., n, a, lda, b, lda, vr, lda,
469 $ alphar, alphai, beta, work,
471 IF( result( 3 ).GT.thresh )
THEN
472 WRITE( nout, fmt = 9998 )
'Right',
'SGGEVX',
473 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
480 IF( s( i ).EQ.zero )
THEN
481 IF( stru( i ).GT.abnorm*ulp )
482 $ result( 3 ) = ulpinv
483 ELSE IF( stru( i ).EQ.zero )
THEN
484 IF( s( i ).GT.abnorm*ulp )
485 $ result( 3 ) = ulpinv
487 work( i ) = max( abs( stru( i ) / s( i ) ),
488 $ abs( s( i ) / stru( i ) ) )
489 result( 3 ) = max( result( 3 ), work( i ) )
496 IF( dif( 1 ).EQ.zero )
THEN
497 IF( diftru( 1 ).GT.abnorm*ulp )
498 $ result( 4 ) = ulpinv
499 ELSE IF( diftru( 1 ).EQ.zero )
THEN
500 IF( dif( 1 ).GT.abnorm*ulp )
501 $ result( 4 ) = ulpinv
502 ELSE IF( dif( 5 ).EQ.zero )
THEN
503 IF( diftru( 5 ).GT.abnorm*ulp )
504 $ result( 4 ) = ulpinv
505 ELSE IF( diftru( 5 ).EQ.zero )
THEN
506 IF( dif( 5 ).GT.abnorm*ulp )
507 $ result( 4 ) = ulpinv
509 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
510 $ abs( dif( 1 ) / diftru( 1 ) ) )
511 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
512 $ abs( dif( 5 ) / diftru( 5 ) ) )
513 result( 4 ) = max( ratio1, ratio2 )
521 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
522 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
528 IF( nerrs.EQ.0 )
THEN
529 WRITE( nout, fmt = 9997 )
'SXV'
535 WRITE( nout, fmt = 9995 )
536 WRITE( nout, fmt = 9994 )
537 WRITE( nout, fmt = 9993 )
541 WRITE( nout, fmt = 9992 )
'''',
546 IF( result( j ).LT.10000.0 )
THEN
547 WRITE( nout, fmt = 9991 )iptype, iwa,
548 $ iwb, iwx, iwy, j, result( j )
550 WRITE( nout, fmt = 9990 )iptype, iwa,
551 $ iwb, iwx, iwy, j, result( j )
571 READ( nin, fmt = *,
END = 150 )n
575 READ( nin, fmt = * )( a( i, j ), j = 1, n )
578 READ( nin, fmt = * )( b( i, j ), j = 1, n )
580 READ( nin, fmt = * )( stru( i ), i = 1, n )
581 READ( nin, fmt = * )( diftru( i ), i = 1, n )
589 CALL
slacpy(
'F', n, n, a, lda, ai, lda )
590 CALL
slacpy(
'F', n, n, b, lda, bi, lda )
592 CALL
sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
593 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
594 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
597 IF( linfo.NE.0 )
THEN
599 WRITE( nout, fmt = 9987 )
'SGGEVX', linfo, n, nptknt
605 CALL
slacpy(
'Full', n, n, ai, lda, work, n )
606 CALL
slacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
607 abnorm =
slange(
'Fro', n, 2*n, work, n, work )
612 CALL
sget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
613 $ beta, work, result( 1 ) )
614 IF( result( 2 ).GT.thresh )
THEN
615 WRITE( nout, fmt = 9986 )
'Left',
'SGGEVX', result( 2 ), n,
620 CALL
sget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
621 $ beta, work, result( 2 ) )
622 IF( result( 3 ).GT.thresh )
THEN
623 WRITE( nout, fmt = 9986 )
'Right',
'SGGEVX', result( 3 ), n,
631 IF( s( i ).EQ.zero )
THEN
632 IF( stru( i ).GT.abnorm*ulp )
633 $ result( 3 ) = ulpinv
634 ELSE IF( stru( i ).EQ.zero )
THEN
635 IF( s( i ).GT.abnorm*ulp )
636 $ result( 3 ) = ulpinv
638 work( i ) = max( abs( stru( i ) / s( i ) ),
639 $ abs( s( i ) / stru( i ) ) )
640 result( 3 ) = max( result( 3 ), work( i ) )
647 IF( dif( 1 ).EQ.zero )
THEN
648 IF( diftru( 1 ).GT.abnorm*ulp )
649 $ result( 4 ) = ulpinv
650 ELSE IF( diftru( 1 ).EQ.zero )
THEN
651 IF( dif( 1 ).GT.abnorm*ulp )
652 $ result( 4 ) = ulpinv
653 ELSE IF( dif( 5 ).EQ.zero )
THEN
654 IF( diftru( 5 ).GT.abnorm*ulp )
655 $ result( 4 ) = ulpinv
656 ELSE IF( diftru( 5 ).EQ.zero )
THEN
657 IF( dif( 5 ).GT.abnorm*ulp )
658 $ result( 4 ) = ulpinv
660 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
661 $ abs( dif( 1 ) / diftru( 1 ) ) )
662 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
663 $ abs( dif( 5 ) / diftru( 5 ) ) )
664 result( 4 ) = max( ratio1, ratio2 )
672 IF( result( j ).GE.thrsh2 )
THEN
677 IF( nerrs.EQ.0 )
THEN
678 WRITE( nout, fmt = 9997 )
'SXV'
684 WRITE( nout, fmt = 9996 )
688 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
692 IF( result( j ).LT.10000.0 )
THEN
693 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
695 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
707 CALL
alasvm(
'SXV', nout, nerrs, ntestt, 0 )
713 9999 format(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
714 $ i6,
', JTYPE=', i6,
')' )
716 9998 format(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
717 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
718 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
719 $
', IWX=', i5,
', IWY=', i5 )
721 9997 format( / 1x, a3,
' -- Real Expert Eigenvalue/vector',
722 $
' problem driver' )
724 9996 format(
' Input Example' )
726 9995 format(
' Matrix types: ', / )
728 9994 format(
' TYPE 1: Da is diagonal, Db is identity, ',
729 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
730 $ /
' YH and X are left and right eigenvectors. ', / )
732 9993 format(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
733 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
734 $ /
' YH and X are left and right eigenvectors. ', / )
736 9992 format( /
' Tests performed: ', / 4x,
737 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
738 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
739 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
740 $ /
' 2 = max | ( b A - a B ) r | / const.',
741 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
742 $
' over all eigenvalues', /
743 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
744 $
' over the 1st and 5th eigenvectors', / )
746 9991 format(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
747 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
748 9990 format(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
749 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, e10.3 )
750 9989 format(
' Input example #', i2,
', matrix order=', i4,
',',
751 $
' result ', i2,
' is', 0p, f8.2 )
752 9988 format(
' Input example #', i2,
', matrix order=', i4,
',',
753 $
' result ', i2,
' is', 1p, e10.3 )
754 9987 format(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
755 $ i6,
', Input example #', i2,
')' )
757 9986 format(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
758 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
759 $
'N=', i6,
', Input Example #', i2,
')' )