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

◆ cget24()

subroutine cget24 ( logical comp,
integer jtype,
real thresh,
integer, dimension( 4 ) iseed,
integer nounit,
integer n,
complex, dimension( lda, * ) a,
integer lda,
complex, dimension( lda, * ) h,
complex, dimension( lda, * ) ht,
complex, dimension( * ) w,
complex, dimension( * ) wt,
complex, dimension( * ) wtmp,
complex, dimension( ldvs, * ) vs,
integer ldvs,
complex, dimension( ldvs, * ) vs1,
real rcdein,
real rcdvin,
integer nslct,
integer, dimension( * ) islct,
integer isrt,
real, dimension( 17 ) result,
complex, dimension( * ) work,
integer lwork,
real, dimension( * ) rwork,
logical, dimension( * ) bwork,
integer info )

CGET24

Purpose:
!>
!>    CGET24 checks the nonsymmetric eigenvalue (Schur form) problem
!>    expert driver CGEESX.
!>
!>    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 W 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 W are eigenvalues of T
!>            1/ulp otherwise
!>            If workspace sufficient, also compare W 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 CGEESX 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 CGEESX 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  if the , 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 COMPLEX 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 COMPLEX array, dimension (LDA, N)
!>          Another copy of the test matrix A, modified by CGEESX.
!>
[out]HT
!>          HT is COMPLEX array, dimension (LDA, N)
!>          Yet another copy of the test matrix A, modified by CGEESX.
!>
[out]W
!>          W is COMPLEX array, dimension (N)
!>          The computed eigenvalues of A.
!>
[out]WT
!>          WT is COMPLEX array, dimension (N)
!>          Like W, this array contains the eigenvalues of A,
!>          but those computed when CGEESX only computes a partial
!>          eigendecomposition, i.e. not Schur vectors
!>
[out]WTMP
!>          WTMP is COMPLEX array, dimension (N)
!>          Like W, this array contains the eigenvalues of A,
!>          but sorted by increasing real or imaginary part.
!>
[out]VS
!>          VS is COMPLEX 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 COMPLEX 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 or imaginary part is
!>          selected. The real part is used if ISRT = 0, and the
!>          imaginary part if ISRT = 1.
!>          Not referenced if COMP = .FALSE.
!>
[in]ISRT
!>          ISRT is INTEGER
!>          When COMP = .TRUE., ISRT describes how ISLCT is used to
!>          choose a subset of the spectrum.
!>          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 COMPLEX array, dimension (2*N*N)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The number of entries in WORK to be passed to CGEESX. This
!>          must be at least 2*N, and N*(N+1)/2 if tests 14--16 are to
!>          be performed.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (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, CGEESX 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 331 of file cget24.f.

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