LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sget24()

subroutine sget24 ( logical  comp,
integer  jtype,
real  thresh,
integer, dimension( 4 )  iseed,
integer  nounit,
integer  n,
real, dimension( lda, * )  a,
integer  lda,
real, dimension( lda, * )  h,
real, dimension( lda, * )  ht,
real, dimension( * )  wr,
real, dimension( * )  wi,
real, dimension( * )  wrt,
real, dimension( * )  wit,
real, dimension( * )  wrtmp,
real, dimension( * )  witmp,
real, dimension( ldvs, * )  vs,
integer  ldvs,
real, dimension( ldvs, * )  vs1,
real  rcdein,
real  rcdvin,
integer  nslct,
integer, dimension( * )  islct,
real, dimension( 17 )  result,
real, dimension( * )  work,
integer  lwork,
integer, dimension( * )  iwork,
logical, dimension( * )  bwork,
integer  info 
)

SGET24

Purpose:
    SGET24 checks the nonsymmetric eigenvalue (Schur form) problem
    expert driver SGEESX.

    If COMP = .FALSE., the first 13 of the following tests will be
    be performed on the input matrix A, and also tests 14 and 15
    if LWORK is sufficiently large.
    If COMP = .TRUE., all 17 test will be performed.

    (1)     0 if T is in Schur form, 1/ulp otherwise
           (no sorting of eigenvalues)

    (2)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (no sorting of eigenvalues).

    (3)     | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).

    (4)     0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            (no sorting of eigenvalues)

    (5)     0     if T(with VS) = T(without VS),
            1/ulp otherwise
            (no sorting of eigenvalues)

    (6)     0     if eigenvalues(with VS) = eigenvalues(without VS),
            1/ulp otherwise
            (no sorting of eigenvalues)

    (7)     0 if T is in Schur form, 1/ulp otherwise
            (with sorting of eigenvalues)

    (8)     | A - VS T VS' | / ( n |A| ulp )

      Here VS is the matrix of Schur eigenvectors, and T is in Schur
      form  (with sorting of eigenvalues).

    (9)     | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).

    (10)    0     if WR+sqrt(-1)*WI are eigenvalues of T
            1/ulp otherwise
            If workspace sufficient, also compare WR, WI with and
            without reciprocal condition numbers
            (with sorting of eigenvalues)

    (11)    0     if T(with VS) = T(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare T with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (12)    0     if eigenvalues(with VS) = eigenvalues(without VS),
            1/ulp otherwise
            If workspace sufficient, also compare VS with and without
            reciprocal condition numbers
            (with sorting of eigenvalues)

    (13)    if sorting worked and SDIM is the number of
            eigenvalues which were SELECTed
            If workspace sufficient, also compare SDIM with and
            without reciprocal condition numbers

    (14)    if RCONDE the same no matter if VS and/or RCONDV computed

    (15)    if RCONDV the same no matter if VS and/or RCONDE computed

    (16)  |RCONDE - RCDEIN| / cond(RCONDE)

       RCONDE is the reciprocal average eigenvalue condition number
       computed by SGEESX and RCDEIN (the precomputed true value)
       is supplied as input.  cond(RCONDE) is the condition number
       of RCONDE, and takes errors in computing RCONDE into account,
       so that the resulting quantity should be O(ULP). cond(RCONDE)
       is essentially given by norm(A)/RCONDV.

    (17)  |RCONDV - RCDVIN| / cond(RCONDV)

       RCONDV is the reciprocal right invariant subspace condition
       number computed by SGEESX and RCDVIN (the precomputed true
       value) is supplied as input. cond(RCONDV) is the condition
       number of RCONDV, and takes errors in computing RCONDV into
       account, so that the resulting quantity should be O(ULP).
       cond(RCONDV) is essentially given by norm(A)/RCONDE.
Parameters
[in]COMP
          COMP is LOGICAL
          COMP describes which input tests to perform:
            = .FALSE. if the computed condition numbers are not to
                      be tested against RCDVIN and RCDEIN
            = .TRUE.  if they are to be compared
[in]JTYPE
          JTYPE is INTEGER
          Type of input matrix. Used to label output if error occurs.
[in]ISEED
          ISEED is INTEGER array, dimension (4)
          If COMP = .FALSE., the random number generator seed
          used to produce matrix.
          If COMP = .TRUE., ISEED(1) = the number of the example.
          Used to label output if error occurs.
[in]THRESH
          THRESH is REAL
          A test will count as "failed" if the "error", computed as
          described above, exceeds THRESH.  Note that the error
          is scaled to be O(1), so THRESH should be a reasonably
          small multiple of 1, e.g., 10 or 100.  In particular,
          it should not depend on the precision (single vs. double)
          or the size of the matrix.  It must be at least zero.
[in]NOUNIT
          NOUNIT is INTEGER
          The FORTRAN unit number for printing out error messages
          (e.g., if a routine returns INFO not equal to 0.)
[in]N
          N is INTEGER
          The dimension of A. N must be at least 0.
[in,out]A
          A is REAL array, dimension (LDA, N)
          Used to hold the matrix whose eigenvalues are to be
          computed.
[in]LDA
          LDA is INTEGER
          The leading dimension of A, and H. LDA must be at
          least 1 and at least N.
[out]H
          H is REAL array, dimension (LDA, N)
          Another copy of the test matrix A, modified by SGEESX.
[out]HT
          HT is REAL array, dimension (LDA, N)
          Yet another copy of the test matrix A, modified by SGEESX.
[out]WR
          WR is REAL array, dimension (N)
[out]WI
          WI is REAL array, dimension (N)

          The real and imaginary parts of the eigenvalues of A.
          On exit, WR + WI*i are the eigenvalues of the matrix in A.
[out]WRT
          WRT is REAL array, dimension (N)
[out]WIT
          WIT is REAL array, dimension (N)

          Like WR, WI, these arrays contain the eigenvalues of A,
          but those computed when SGEESX only computes a partial
          eigendecomposition, i.e. not Schur vectors
[out]WRTMP
          WRTMP is REAL array, dimension (N)
[out]WITMP
          WITMP is REAL array, dimension (N)

          Like WR, WI, these arrays contain the eigenvalues of A,
          but sorted by increasing real part.
[out]VS
          VS is REAL array, dimension (LDVS, N)
          VS holds the computed Schur vectors.
[in]LDVS
          LDVS is INTEGER
          Leading dimension of VS. Must be at least max(1, N).
[out]VS1
          VS1 is REAL array, dimension (LDVS, N)
          VS1 holds another copy of the computed Schur vectors.
[in]RCDEIN
          RCDEIN is REAL
          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
          condition number for the average of selected eigenvalues.
[in]RCDVIN
          RCDVIN is REAL
          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
          condition number for the selected right invariant subspace.
[in]NSLCT
          NSLCT is INTEGER
          When COMP = .TRUE. the number of selected eigenvalues
          corresponding to the precomputed values RCDEIN and RCDVIN.
[in]ISLCT
          ISLCT is INTEGER array, dimension (NSLCT)
          When COMP = .TRUE. ISLCT selects the eigenvalues of the
          input matrix corresponding to the precomputed values RCDEIN
          and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
          eigenvalue with the J-th largest real part is selected.
          Not referenced if COMP = .FALSE.
[out]RESULT
          RESULT is REAL array, dimension (17)
          The values computed by the 17 tests described above.
          The values are currently limited to 1/ulp, to avoid
          overflow.
[out]WORK
          WORK is REAL array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The number of entries in WORK to be passed to SGEESX. This
          must be at least 3*N, and N+N**2 if tests 14--16 are to
          be performed.
[out]IWORK
          IWORK is INTEGER array, dimension (N*N)
[out]BWORK
          BWORK is LOGICAL array, dimension (N)
[out]INFO
          INFO is INTEGER
          If 0,  successful exit.
          If <0, input parameter -INFO had an incorrect value.
          If >0, SGEESX returned an error code, the absolute
                 value of which is returned.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 339 of file sget24.f.

343*
344* -- LAPACK test routine --
345* -- LAPACK is a software package provided by Univ. of Tennessee, --
346* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
347*
348* .. Scalar Arguments ..
349 LOGICAL COMP
350 INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
351 REAL RCDEIN, RCDVIN, THRESH
352* ..
353* .. Array Arguments ..
354 LOGICAL BWORK( * )
355 INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
356 REAL A( LDA, * ), H( LDA, * ), HT( LDA, * ),
357 $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
358 $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
359 $ WR( * ), WRT( * ), WRTMP( * )
360* ..
361*
362* =====================================================================
363*
364* .. Parameters ..
365 REAL ZERO, ONE
366 parameter( zero = 0.0e0, one = 1.0e0 )
367 REAL EPSIN
368 parameter( epsin = 5.9605e-8 )
369* ..
370* .. Local Scalars ..
371 CHARACTER SORT
372 INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
373 $ RSUB, SDIM, SDIM1
374 REAL ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
375 $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
376 $ VRMIN, WNORM
377* ..
378* .. Local Arrays ..
379 INTEGER IPNT( 20 )
380* ..
381* .. Arrays in Common ..
382 LOGICAL SELVAL( 20 )
383 REAL SELWI( 20 ), SELWR( 20 )
384* ..
385* .. Scalars in Common ..
386 INTEGER SELDIM, SELOPT
387* ..
388* .. Common blocks ..
389 COMMON / sslct / selopt, seldim, selval, selwr, selwi
390* ..
391* .. External Functions ..
392 LOGICAL SSLECT
393 REAL SLAMCH, SLANGE
394 EXTERNAL sslect, slamch, slange
395* ..
396* .. External Subroutines ..
397 EXTERNAL scopy, sgeesx, sgemm, slacpy, sort01, xerbla
398* ..
399* .. Intrinsic Functions ..
400 INTRINSIC abs, max, min, real, sign, sqrt
401* ..
402* .. Executable Statements ..
403*
404* Check for errors
405*
406 info = 0
407 IF( thresh.LT.zero ) THEN
408 info = -3
409 ELSE IF( nounit.LE.0 ) THEN
410 info = -5
411 ELSE IF( n.LT.0 ) THEN
412 info = -6
413 ELSE IF( lda.LT.1 .OR. lda.LT.n ) THEN
414 info = -8
415 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.n ) THEN
416 info = -18
417 ELSE IF( lwork.LT.3*n ) THEN
418 info = -26
419 END IF
420*
421 IF( info.NE.0 ) THEN
422 CALL xerbla( 'SGET24', -info )
423 RETURN
424 END IF
425*
426* Quick return if nothing to do
427*
428 DO 10 i = 1, 17
429 result( i ) = -one
430 10 CONTINUE
431*
432 IF( n.EQ.0 )
433 $ RETURN
434*
435* Important constants
436*
437 smlnum = slamch( 'Safe minimum' )
438 ulp = slamch( 'Precision' )
439 ulpinv = one / ulp
440*
441* Perform tests (1)-(13)
442*
443 selopt = 0
444 liwork = n*n
445 DO 120 isort = 0, 1
446 IF( isort.EQ.0 ) THEN
447 sort = 'N'
448 rsub = 0
449 ELSE
450 sort = 'S'
451 rsub = 6
452 END IF
453*
454* Compute Schur form and Schur vectors, and test them
455*
456 CALL slacpy( 'F', n, n, a, lda, h, lda )
457 CALL sgeesx( 'V', sort, sslect, 'N', n, h, lda, sdim, wr, wi,
458 $ vs, ldvs, rconde, rcondv, work, lwork, iwork,
459 $ liwork, bwork, iinfo )
460 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
461 result( 1+rsub ) = ulpinv
462 IF( jtype.NE.22 ) THEN
463 WRITE( nounit, fmt = 9998 )'SGEESX1', iinfo, n, jtype,
464 $ iseed
465 ELSE
466 WRITE( nounit, fmt = 9999 )'SGEESX1', iinfo, n,
467 $ iseed( 1 )
468 END IF
469 info = abs( iinfo )
470 RETURN
471 END IF
472 IF( isort.EQ.0 ) THEN
473 CALL scopy( n, wr, 1, wrtmp, 1 )
474 CALL scopy( n, wi, 1, witmp, 1 )
475 END IF
476*
477* Do Test (1) or Test (7)
478*
479 result( 1+rsub ) = zero
480 DO 30 j = 1, n - 2
481 DO 20 i = j + 2, n
482 IF( h( i, j ).NE.zero )
483 $ result( 1+rsub ) = ulpinv
484 20 CONTINUE
485 30 CONTINUE
486 DO 40 i = 1, n - 2
487 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.zero )
488 $ result( 1+rsub ) = ulpinv
489 40 CONTINUE
490 DO 50 i = 1, n - 1
491 IF( h( i+1, i ).NE.zero ) THEN
492 IF( h( i, i ).NE.h( i+1, i+1 ) .OR. h( i, i+1 ).EQ.
493 $ zero .OR. sign( one, h( i+1, i ) ).EQ.
494 $ sign( one, h( i, i+1 ) ) )result( 1+rsub ) = ulpinv
495 END IF
496 50 CONTINUE
497*
498* Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
499*
500* Copy A to VS1, used as workspace
501*
502 CALL slacpy( ' ', n, n, a, lda, vs1, ldvs )
503*
504* Compute Q*H and store in HT.
505*
506 CALL sgemm( 'No transpose', 'No transpose', n, n, n, one, vs,
507 $ ldvs, h, lda, zero, ht, lda )
508*
509* Compute A - Q*H*Q'
510*
511 CALL sgemm( 'No transpose', 'Transpose', n, n, n, -one, ht,
512 $ lda, vs, ldvs, one, vs1, ldvs )
513*
514 anorm = max( slange( '1', n, n, a, lda, work ), smlnum )
515 wnorm = slange( '1', n, n, vs1, ldvs, work )
516*
517 IF( anorm.GT.wnorm ) THEN
518 result( 2+rsub ) = ( wnorm / anorm ) / ( n*ulp )
519 ELSE
520 IF( anorm.LT.one ) THEN
521 result( 2+rsub ) = ( min( wnorm, n*anorm ) / anorm ) /
522 $ ( n*ulp )
523 ELSE
524 result( 2+rsub ) = min( wnorm / anorm, real( n ) ) /
525 $ ( n*ulp )
526 END IF
527 END IF
528*
529* Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
530*
531 CALL sort01( 'Columns', n, n, vs, ldvs, work, lwork,
532 $ result( 3+rsub ) )
533*
534* Do Test (4) or Test (10)
535*
536 result( 4+rsub ) = zero
537 DO 60 i = 1, n
538 IF( h( i, i ).NE.wr( i ) )
539 $ result( 4+rsub ) = ulpinv
540 60 CONTINUE
541 IF( n.GT.1 ) THEN
542 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
543 $ result( 4+rsub ) = ulpinv
544 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
545 $ result( 4+rsub ) = ulpinv
546 END IF
547 DO 70 i = 1, n - 1
548 IF( h( i+1, i ).NE.zero ) THEN
549 tmp = sqrt( abs( h( i+1, i ) ) )*
550 $ sqrt( abs( h( i, i+1 ) ) )
551 result( 4+rsub ) = max( result( 4+rsub ),
552 $ abs( wi( i )-tmp ) /
553 $ max( ulp*tmp, smlnum ) )
554 result( 4+rsub ) = max( result( 4+rsub ),
555 $ abs( wi( i+1 )+tmp ) /
556 $ max( ulp*tmp, smlnum ) )
557 ELSE IF( i.GT.1 ) THEN
558 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.zero .AND.
559 $ wi( i ).NE.zero )result( 4+rsub ) = ulpinv
560 END IF
561 70 CONTINUE
562*
563* Do Test (5) or Test (11)
564*
565 CALL slacpy( 'F', n, n, a, lda, ht, lda )
566 CALL sgeesx( 'N', sort, sslect, 'N', n, ht, lda, sdim, wrt,
567 $ wit, vs, ldvs, rconde, rcondv, work, lwork, iwork,
568 $ liwork, bwork, iinfo )
569 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
570 result( 5+rsub ) = ulpinv
571 IF( jtype.NE.22 ) THEN
572 WRITE( nounit, fmt = 9998 )'SGEESX2', iinfo, n, jtype,
573 $ iseed
574 ELSE
575 WRITE( nounit, fmt = 9999 )'SGEESX2', iinfo, n,
576 $ iseed( 1 )
577 END IF
578 info = abs( iinfo )
579 GO TO 250
580 END IF
581*
582 result( 5+rsub ) = zero
583 DO 90 j = 1, n
584 DO 80 i = 1, n
585 IF( h( i, j ).NE.ht( i, j ) )
586 $ result( 5+rsub ) = ulpinv
587 80 CONTINUE
588 90 CONTINUE
589*
590* Do Test (6) or Test (12)
591*
592 result( 6+rsub ) = zero
593 DO 100 i = 1, n
594 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
595 $ result( 6+rsub ) = ulpinv
596 100 CONTINUE
597*
598* Do Test (13)
599*
600 IF( isort.EQ.1 ) THEN
601 result( 13 ) = zero
602 knteig = 0
603 DO 110 i = 1, n
604 IF( sslect( wr( i ), wi( i ) ) .OR.
605 $ sslect( wr( i ), -wi( i ) ) )knteig = knteig + 1
606 IF( i.LT.n ) THEN
607 IF( ( sslect( wr( i+1 ), wi( i+1 ) ) .OR.
608 $ sslect( wr( i+1 ), -wi( i+1 ) ) ) .AND.
609 $ ( .NOT.( sslect( wr( i ),
610 $ wi( i ) ) .OR. sslect( wr( i ),
611 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )result( 13 )
612 $ = ulpinv
613 END IF
614 110 CONTINUE
615 IF( sdim.NE.knteig )
616 $ result( 13 ) = ulpinv
617 END IF
618*
619 120 CONTINUE
620*
621* If there is enough workspace, perform tests (14) and (15)
622* as well as (10) through (13)
623*
624 IF( lwork.GE.n+( n*n ) / 2 ) THEN
625*
626* Compute both RCONDE and RCONDV with VS
627*
628 sort = 'S'
629 result( 14 ) = zero
630 result( 15 ) = zero
631 CALL slacpy( 'F', n, n, a, lda, ht, lda )
632 CALL sgeesx( 'V', sort, sslect, 'B', n, ht, lda, sdim1, wrt,
633 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
634 $ iwork, liwork, bwork, iinfo )
635 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
636 result( 14 ) = ulpinv
637 result( 15 ) = ulpinv
638 IF( jtype.NE.22 ) THEN
639 WRITE( nounit, fmt = 9998 )'SGEESX3', iinfo, n, jtype,
640 $ iseed
641 ELSE
642 WRITE( nounit, fmt = 9999 )'SGEESX3', iinfo, n,
643 $ iseed( 1 )
644 END IF
645 info = abs( iinfo )
646 GO TO 250
647 END IF
648*
649* Perform tests (10), (11), (12), and (13)
650*
651 DO 140 i = 1, n
652 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
653 $ result( 10 ) = ulpinv
654 DO 130 j = 1, n
655 IF( h( i, j ).NE.ht( i, j ) )
656 $ result( 11 ) = ulpinv
657 IF( vs( i, j ).NE.vs1( i, j ) )
658 $ result( 12 ) = ulpinv
659 130 CONTINUE
660 140 CONTINUE
661 IF( sdim.NE.sdim1 )
662 $ result( 13 ) = ulpinv
663*
664* Compute both RCONDE and RCONDV without VS, and compare
665*
666 CALL slacpy( 'F', n, n, a, lda, ht, lda )
667 CALL sgeesx( 'N', sort, sslect, 'B', n, ht, lda, sdim1, wrt,
668 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
669 $ iwork, liwork, bwork, iinfo )
670 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
671 result( 14 ) = ulpinv
672 result( 15 ) = ulpinv
673 IF( jtype.NE.22 ) THEN
674 WRITE( nounit, fmt = 9998 )'SGEESX4', iinfo, n, jtype,
675 $ iseed
676 ELSE
677 WRITE( nounit, fmt = 9999 )'SGEESX4', iinfo, n,
678 $ iseed( 1 )
679 END IF
680 info = abs( iinfo )
681 GO TO 250
682 END IF
683*
684* Perform tests (14) and (15)
685*
686 IF( rcnde1.NE.rconde )
687 $ result( 14 ) = ulpinv
688 IF( rcndv1.NE.rcondv )
689 $ result( 15 ) = ulpinv
690*
691* Perform tests (10), (11), (12), and (13)
692*
693 DO 160 i = 1, n
694 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
695 $ result( 10 ) = ulpinv
696 DO 150 j = 1, n
697 IF( h( i, j ).NE.ht( i, j ) )
698 $ result( 11 ) = ulpinv
699 IF( vs( i, j ).NE.vs1( i, j ) )
700 $ result( 12 ) = ulpinv
701 150 CONTINUE
702 160 CONTINUE
703 IF( sdim.NE.sdim1 )
704 $ result( 13 ) = ulpinv
705*
706* Compute RCONDE with VS, and compare
707*
708 CALL slacpy( 'F', n, n, a, lda, ht, lda )
709 CALL sgeesx( 'V', sort, sslect, 'E', n, ht, lda, sdim1, wrt,
710 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
711 $ iwork, liwork, bwork, iinfo )
712 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
713 result( 14 ) = ulpinv
714 IF( jtype.NE.22 ) THEN
715 WRITE( nounit, fmt = 9998 )'SGEESX5', iinfo, n, jtype,
716 $ iseed
717 ELSE
718 WRITE( nounit, fmt = 9999 )'SGEESX5', iinfo, n,
719 $ iseed( 1 )
720 END IF
721 info = abs( iinfo )
722 GO TO 250
723 END IF
724*
725* Perform test (14)
726*
727 IF( rcnde1.NE.rconde )
728 $ result( 14 ) = ulpinv
729*
730* Perform tests (10), (11), (12), and (13)
731*
732 DO 180 i = 1, n
733 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
734 $ result( 10 ) = ulpinv
735 DO 170 j = 1, n
736 IF( h( i, j ).NE.ht( i, j ) )
737 $ result( 11 ) = ulpinv
738 IF( vs( i, j ).NE.vs1( i, j ) )
739 $ result( 12 ) = ulpinv
740 170 CONTINUE
741 180 CONTINUE
742 IF( sdim.NE.sdim1 )
743 $ result( 13 ) = ulpinv
744*
745* Compute RCONDE without VS, and compare
746*
747 CALL slacpy( 'F', n, n, a, lda, ht, lda )
748 CALL sgeesx( 'N', sort, sslect, 'E', n, ht, lda, sdim1, wrt,
749 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
750 $ iwork, liwork, bwork, iinfo )
751 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
752 result( 14 ) = ulpinv
753 IF( jtype.NE.22 ) THEN
754 WRITE( nounit, fmt = 9998 )'SGEESX6', iinfo, n, jtype,
755 $ iseed
756 ELSE
757 WRITE( nounit, fmt = 9999 )'SGEESX6', iinfo, n,
758 $ iseed( 1 )
759 END IF
760 info = abs( iinfo )
761 GO TO 250
762 END IF
763*
764* Perform test (14)
765*
766 IF( rcnde1.NE.rconde )
767 $ result( 14 ) = ulpinv
768*
769* Perform tests (10), (11), (12), and (13)
770*
771 DO 200 i = 1, n
772 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
773 $ result( 10 ) = ulpinv
774 DO 190 j = 1, n
775 IF( h( i, j ).NE.ht( i, j ) )
776 $ result( 11 ) = ulpinv
777 IF( vs( i, j ).NE.vs1( i, j ) )
778 $ result( 12 ) = ulpinv
779 190 CONTINUE
780 200 CONTINUE
781 IF( sdim.NE.sdim1 )
782 $ result( 13 ) = ulpinv
783*
784* Compute RCONDV with VS, and compare
785*
786 CALL slacpy( 'F', n, n, a, lda, ht, lda )
787 CALL sgeesx( 'V', sort, sslect, 'V', n, ht, lda, sdim1, wrt,
788 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
789 $ iwork, liwork, bwork, iinfo )
790 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
791 result( 15 ) = ulpinv
792 IF( jtype.NE.22 ) THEN
793 WRITE( nounit, fmt = 9998 )'SGEESX7', iinfo, n, jtype,
794 $ iseed
795 ELSE
796 WRITE( nounit, fmt = 9999 )'SGEESX7', iinfo, n,
797 $ iseed( 1 )
798 END IF
799 info = abs( iinfo )
800 GO TO 250
801 END IF
802*
803* Perform test (15)
804*
805 IF( rcndv1.NE.rcondv )
806 $ result( 15 ) = ulpinv
807*
808* Perform tests (10), (11), (12), and (13)
809*
810 DO 220 i = 1, n
811 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
812 $ result( 10 ) = ulpinv
813 DO 210 j = 1, n
814 IF( h( i, j ).NE.ht( i, j ) )
815 $ result( 11 ) = ulpinv
816 IF( vs( i, j ).NE.vs1( i, j ) )
817 $ result( 12 ) = ulpinv
818 210 CONTINUE
819 220 CONTINUE
820 IF( sdim.NE.sdim1 )
821 $ result( 13 ) = ulpinv
822*
823* Compute RCONDV without VS, and compare
824*
825 CALL slacpy( 'F', n, n, a, lda, ht, lda )
826 CALL sgeesx( 'N', sort, sslect, 'V', n, ht, lda, sdim1, wrt,
827 $ wit, vs1, ldvs, rcnde1, rcndv1, work, lwork,
828 $ iwork, liwork, bwork, iinfo )
829 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
830 result( 15 ) = ulpinv
831 IF( jtype.NE.22 ) THEN
832 WRITE( nounit, fmt = 9998 )'SGEESX8', iinfo, n, jtype,
833 $ iseed
834 ELSE
835 WRITE( nounit, fmt = 9999 )'SGEESX8', iinfo, n,
836 $ iseed( 1 )
837 END IF
838 info = abs( iinfo )
839 GO TO 250
840 END IF
841*
842* Perform test (15)
843*
844 IF( rcndv1.NE.rcondv )
845 $ result( 15 ) = ulpinv
846*
847* Perform tests (10), (11), (12), and (13)
848*
849 DO 240 i = 1, n
850 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
851 $ result( 10 ) = ulpinv
852 DO 230 j = 1, n
853 IF( h( i, j ).NE.ht( i, j ) )
854 $ result( 11 ) = ulpinv
855 IF( vs( i, j ).NE.vs1( i, j ) )
856 $ result( 12 ) = ulpinv
857 230 CONTINUE
858 240 CONTINUE
859 IF( sdim.NE.sdim1 )
860 $ result( 13 ) = ulpinv
861*
862 END IF
863*
864 250 CONTINUE
865*
866* If there are precomputed reciprocal condition numbers, compare
867* computed values with them.
868*
869 IF( comp ) THEN
870*
871* First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
872* the logical function SSLECT selects the eigenvalues specified
873* by NSLCT and ISLCT.
874*
875 seldim = n
876 selopt = 1
877 eps = max( ulp, epsin )
878 DO 260 i = 1, n
879 ipnt( i ) = i
880 selval( i ) = .false.
881 selwr( i ) = wrtmp( i )
882 selwi( i ) = witmp( i )
883 260 CONTINUE
884 DO 280 i = 1, n - 1
885 kmin = i
886 vrmin = wrtmp( i )
887 vimin = witmp( i )
888 DO 270 j = i + 1, n
889 IF( wrtmp( j ).LT.vrmin ) THEN
890 kmin = j
891 vrmin = wrtmp( j )
892 vimin = witmp( j )
893 END IF
894 270 CONTINUE
895 wrtmp( kmin ) = wrtmp( i )
896 witmp( kmin ) = witmp( i )
897 wrtmp( i ) = vrmin
898 witmp( i ) = vimin
899 itmp = ipnt( i )
900 ipnt( i ) = ipnt( kmin )
901 ipnt( kmin ) = itmp
902 280 CONTINUE
903 DO 290 i = 1, nslct
904 selval( ipnt( islct( i ) ) ) = .true.
905 290 CONTINUE
906*
907* Compute condition numbers
908*
909 CALL slacpy( 'F', n, n, a, lda, ht, lda )
910 CALL sgeesx( 'N', 'S', sslect, 'B', n, ht, lda, sdim1, wrt,
911 $ wit, vs1, ldvs, rconde, rcondv, work, lwork,
912 $ iwork, liwork, bwork, iinfo )
913 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 ) THEN
914 result( 16 ) = ulpinv
915 result( 17 ) = ulpinv
916 WRITE( nounit, fmt = 9999 )'SGEESX9', iinfo, n, iseed( 1 )
917 info = abs( iinfo )
918 GO TO 300
919 END IF
920*
921* Compare condition number for average of selected eigenvalues
922* taking its condition number into account
923*
924 anorm = slange( '1', n, n, a, lda, work )
925 v = max( real( n )*eps*anorm, smlnum )
926 IF( anorm.EQ.zero )
927 $ v = one
928 IF( v.GT.rcondv ) THEN
929 tol = one
930 ELSE
931 tol = v / rcondv
932 END IF
933 IF( v.GT.rcdvin ) THEN
934 tolin = one
935 ELSE
936 tolin = v / rcdvin
937 END IF
938 tol = max( tol, smlnum / eps )
939 tolin = max( tolin, smlnum / eps )
940 IF( eps*( rcdein-tolin ).GT.rconde+tol ) THEN
941 result( 16 ) = ulpinv
942 ELSE IF( rcdein-tolin.GT.rconde+tol ) THEN
943 result( 16 ) = ( rcdein-tolin ) / ( rconde+tol )
944 ELSE IF( rcdein+tolin.LT.eps*( rconde-tol ) ) THEN
945 result( 16 ) = ulpinv
946 ELSE IF( rcdein+tolin.LT.rconde-tol ) THEN
947 result( 16 ) = ( rconde-tol ) / ( rcdein+tolin )
948 ELSE
949 result( 16 ) = one
950 END IF
951*
952* Compare condition numbers for right invariant subspace
953* taking its condition number into account
954*
955 IF( v.GT.rcondv*rconde ) THEN
956 tol = rcondv
957 ELSE
958 tol = v / rconde
959 END IF
960 IF( v.GT.rcdvin*rcdein ) THEN
961 tolin = rcdvin
962 ELSE
963 tolin = v / rcdein
964 END IF
965 tol = max( tol, smlnum / eps )
966 tolin = max( tolin, smlnum / eps )
967 IF( eps*( rcdvin-tolin ).GT.rcondv+tol ) THEN
968 result( 17 ) = ulpinv
969 ELSE IF( rcdvin-tolin.GT.rcondv+tol ) THEN
970 result( 17 ) = ( rcdvin-tolin ) / ( rcondv+tol )
971 ELSE IF( rcdvin+tolin.LT.eps*( rcondv-tol ) ) THEN
972 result( 17 ) = ulpinv
973 ELSE IF( rcdvin+tolin.LT.rcondv-tol ) THEN
974 result( 17 ) = ( rcondv-tol ) / ( rcdvin+tolin )
975 ELSE
976 result( 17 ) = one
977 END IF
978*
979 300 CONTINUE
980*
981 END IF
982*
983 9999 FORMAT( ' SGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
984 $ i6, ', INPUT EXAMPLE NUMBER = ', i4 )
985 9998 FORMAT( ' SGET24: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
986 $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
987*
988 RETURN
989*
990* End of SGET24
991*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgeesx(jobvs, sort, select, sense, n, a, lda, sdim, wr, wi, vs, ldvs, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
SGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition sgeesx.f:281
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:114
subroutine sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01
Definition sort01.f:116
logical function sslect(zr, zi)
SSLECT
Definition sslect.f:62
Here is the call graph for this function:
Here is the caller graph for this function: