sget23.f

来自「famous linear algebra library (LAPACK) p」· F 代码 · 共 722 行 · 第 1/2 页

      SUBROUTINE SGET23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
     $                   A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
     $                   LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
     $                   RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
     $                   WORK, LWORK, IWORK, INFO )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      LOGICAL            COMP
      CHARACTER          BALANC
      INTEGER            INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
     $                   NOUNIT
      REAL               THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            ISEED( 4 ), IWORK( * )
      REAL               A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
     $                   RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
     $                   RCNDV1( * ), RCONDE( * ), RCONDV( * ),
     $                   RESULT( 11 ), SCALE( * ), SCALE1( * ),
     $                   VL( LDVL, * ), VR( LDVR, * ), WI( * ),
     $                   WI1( * ), WORK( * ), WR( * ), WR1( * )
*     ..
*
*  Purpose
*  =======
*
*     SGET23  checks the nonsymmetric eigenvalue problem driver SGEEVX.
*     If COMP = .FALSE., the first 8 of the following tests will be
*     performed on the input matrix A, and also test 9 if LWORK is
*     sufficiently large.
*     if COMP is .TRUE. all 11 tests will be performed.
*
*     (1)     | A * VR - VR * W | / ( n |A| ulp )
*
*       Here VR is the matrix of unit right eigenvectors.
*       W is a block diagonal matrix, with a 1x1 block for each
*       real eigenvalue and a 2x2 block for each complex conjugate
*       pair.  If eigenvalues j and j+1 are a complex conjugate pair,
*       so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
*       2 x 2 block corresponding to the pair will be:
*
*               (  wr  wi  )
*               ( -wi  wr  )
*
*       Such a block multiplying an n x 2 matrix  ( ur ui ) on the
*       right will be the same as multiplying  ur + i*ui  by  wr + i*wi.
*
*     (2)     | A**H * VL - VL * W**H | / ( n |A| ulp )
*
*       Here VL is the matrix of unit left eigenvectors, A**H is the
*       conjugate transpose of A, and W is as above.
*
*     (3)     | |VR(i)| - 1 | / ulp and largest component real
*
*       VR(i) denotes the i-th column of VR.
*
*     (4)     | |VL(i)| - 1 | / ulp and largest component real
*
*       VL(i) denotes the i-th column of VL.
*
*     (5)     0 if W(full) = W(partial), 1/ulp otherwise
*
*       W(full) denotes the eigenvalues computed when VR, VL, RCONDV
*       and RCONDE are also computed, and W(partial) denotes the
*       eigenvalues computed when only some of VR, VL, RCONDV, and
*       RCONDE are computed.
*
*     (6)     0 if VR(full) = VR(partial), 1/ulp otherwise
*
*       VR(full) denotes the right eigenvectors computed when VL, RCONDV
*       and RCONDE are computed, and VR(partial) denotes the result
*       when only some of VL and RCONDV are computed.
*
*     (7)     0 if VL(full) = VL(partial), 1/ulp otherwise
*
*       VL(full) denotes the left eigenvectors computed when VR, RCONDV
*       and RCONDE are computed, and VL(partial) denotes the result
*       when only some of VR and RCONDV are computed.
*
*     (8)     0 if SCALE, ILO, IHI, ABNRM (full) =
*                  SCALE, ILO, IHI, ABNRM (partial)
*             1/ulp otherwise
*
*       SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
*       (full) is when VR, VL, RCONDE and RCONDV are also computed, and
*       (partial) is when some are not computed.
*
*     (9)     0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise
*
*       RCONDV(full) denotes the reciprocal condition numbers of the
*       right eigenvectors computed when VR, VL and RCONDE are also
*       computed. RCONDV(partial) denotes the reciprocal condition
*       numbers when only some of VR, VL and RCONDE are computed.
*
*    (10)     |RCONDV - RCDVIN| / cond(RCONDV)
*
*       RCONDV is the reciprocal right eigenvector condition number
*       computed by SGEEVX 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.
*
*    (11)     |RCONDE - RCDEIN| / cond(RCONDE)
*
*       RCONDE is the reciprocal eigenvalue condition number
*       computed by SGEEVX 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.
*
*  Arguments
*  =========
*
*  COMP    (input) 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
*
*  BALANC  (input) CHARACTER
*          Describes the balancing option to be tested.
*            = 'N' for no permuting or diagonal scaling
*            = 'P' for permuting but no diagonal scaling
*            = 'S' for no permuting but diagonal scaling
*            = 'B' for permuting and diagonal scaling
*
*  JTYPE   (input) INTEGER
*          Type of input matrix. Used to label output if error occurs.
*
*  THRESH  (input) 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.
*
*  ISEED   (input) 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.
*
*  NOUNIT  (input) INTEGER
*          The FORTRAN unit number for printing out error messages
*          (e.g., if a routine returns INFO not equal to 0.)
*
*  N       (input) INTEGER
*          The dimension of A. N must be at least 0.
*
*  A       (input/output) REAL array, dimension (LDA,N)
*          Used to hold the matrix whose eigenvalues are to be
*          computed.
*
*  LDA     (input) INTEGER
*          The leading dimension of A, and H. LDA must be at
*          least 1 and at least N.
*
*  H       (workspace) REAL array, dimension (LDA,N)
*          Another copy of the test matrix A, modified by SGEEVX.
*
*  WR      (workspace) REAL array, dimension (N)
*  WI      (workspace) 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.
*
*  WR1     (workspace) REAL array, dimension (N)
*  WI1     (workspace) REAL array, dimension (N)
*          Like WR, WI, these arrays contain the eigenvalues of A,
*          but those computed when SGEEVX only computes a partial
*          eigendecomposition, i.e. not the eigenvalues and left
*          and right eigenvectors.
*
*  VL      (workspace) REAL array, dimension (LDVL,N)
*          VL holds the computed left eigenvectors.
*
*  LDVL    (input) INTEGER
*          Leading dimension of VL. Must be at least max(1,N).
*
*  VR      (workspace) REAL array, dimension (LDVR,N)
*          VR holds the computed right eigenvectors.
*
*  LDVR    (input) INTEGER
*          Leading dimension of VR. Must be at least max(1,N).
*
*  LRE     (workspace) REAL array, dimension (LDLRE,N)
*          LRE holds the computed right or left eigenvectors.
*
*  LDLRE   (input) INTEGER
*          Leading dimension of LRE. Must be at least max(1,N).
*
*  RCONDV  (workspace) REAL array, dimension (N)
*          RCONDV holds the computed reciprocal condition numbers
*          for eigenvectors.
*
*  RCNDV1  (workspace) REAL array, dimension (N)
*          RCNDV1 holds more computed reciprocal condition numbers
*          for eigenvectors.
*
*  RCDVIN  (input) REAL array, dimension (N)
*          When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
*          condition numbers for eigenvectors to be compared with
*          RCONDV.
*
*  RCONDE  (workspace) REAL array, dimension (N)
*          RCONDE holds the computed reciprocal condition numbers
*          for eigenvalues.
*
*  RCNDE1  (workspace) REAL array, dimension (N)
*          RCNDE1 holds more computed reciprocal condition numbers
*          for eigenvalues.
*
*  RCDEIN  (input) REAL array, dimension (N)
*          When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
*          condition numbers for eigenvalues to be compared with
*          RCONDE.
*
*  SCALE   (workspace) REAL array, dimension (N)
*          Holds information describing balancing of matrix.
*
*  SCALE1  (workspace) REAL array, dimension (N)
*          Holds information describing balancing of matrix.
*
*  RESULT  (output) REAL array, dimension (11)
*          The values computed by the 11 tests described above.
*          The values are currently limited to 1/ulp, to avoid
*          overflow.
*
*  WORK    (workspace) REAL array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The number of entries in WORK.  This must be at least
*          3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed.
*
*  IWORK   (workspace) INTEGER array, dimension (2*N)
*
*  INFO    (output) INTEGER
*          If 0,  successful exit.
*          If <0, input parameter -INFO had an incorrect value.
*          If >0, SGEEVX returned an error code, the absolute
*                 value of which is returned.
*
*  =====================================================================
*
*
*     .. Parameters ..
      REAL               ZERO, ONE, TWO
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
      REAL               EPSIN
      PARAMETER          ( EPSIN = 5.9605E-8 )
*     ..
*     .. Local Scalars ..
      LOGICAL            BALOK, NOBAL
      CHARACTER          SENSE
      INTEGER            I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
     $                   J, JJ, KMIN
      REAL               ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
     $                   ULP, ULPINV, V, VIMIN, VMAX, VMX, VRMIN, VRMX,
     $                   VTST
*     ..
*     .. Local Arrays ..
      CHARACTER          SENS( 2 )
      REAL               DUM( 1 ), RES( 2 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH, SLAPY2, SNRM2
      EXTERNAL           LSAME, SLAMCH, SLAPY2, SNRM2
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEEVX, SGET22, SLACPY, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN, REAL
*     ..
*     .. Data statements ..
      DATA               SENS / 'N', 'V' /
*     ..
*     .. Executable Statements ..
*
*     Check for errors
*
      NOBAL = LSAME( BALANC, 'N' )
      BALOK = NOBAL .OR. LSAME( BALANC, 'P' ) .OR.
     $        LSAME( BALANC, 'S' ) .OR. LSAME( BALANC, 'B' )
      INFO = 0
      IF( .NOT.BALOK ) THEN
         INFO = -2
      ELSE IF( THRESH.LT.ZERO ) THEN
         INFO = -4
      ELSE IF( NOUNIT.LE.0 ) THEN
         INFO = -6
      ELSE IF( N.LT.0 ) THEN
         INFO = -7
      ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
         INFO = -9
      ELSE IF( LDVL.LT.1 .OR. LDVL.LT.N ) THEN
         INFO = -16
      ELSE IF( LDVR.LT.1 .OR. LDVR.LT.N ) THEN
         INFO = -18
      ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.N ) THEN
         INFO = -20
      ELSE IF( LWORK.LT.3*N .OR. ( COMP .AND. LWORK.LT.6*N+N*N ) ) THEN
         INFO = -31
      END IF
*
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SGET23', -INFO )
         RETURN
      END IF
*
*     Quick return if nothing to do
*
      DO 10 I = 1, 11
         RESULT( I ) = -ONE
   10 CONTINUE
*
      IF( N.EQ.0 )
     $   RETURN
*
*     More Important constants
*
      ULP = SLAMCH( 'Precision' )
      SMLNUM = SLAMCH( 'S' )
      ULPINV = ONE / ULP
*
*     Compute eigenvalues and eigenvectors, and test them
*
      IF( LWORK.GE.6*N+N*N ) THEN
         SENSE = 'B'
         ISENSM = 2
      ELSE
         SENSE = 'E'
         ISENSM = 1
      END IF
      CALL SLACPY( 'F', N, N, A, LDA, H, LDA )
      CALL SGEEVX( BALANC, 'V', 'V', SENSE, N, H, LDA, WR, WI, VL, LDVL,
     $             VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
     $             WORK, LWORK, IWORK, IINFO )
      IF( IINFO.NE.0 ) THEN
         RESULT( 1 ) = ULPINV
         IF( JTYPE.NE.22 ) THEN
            WRITE( NOUNIT, FMT = 9998 )'SGEEVX1', IINFO, N, JTYPE,
     $         BALANC, ISEED
         ELSE
            WRITE( NOUNIT, FMT = 9999 )'SGEEVX1', IINFO, N, ISEED( 1 )
         END IF
         INFO = ABS( IINFO )
         RETURN
      END IF
*
*     Do Test (1)
*
      CALL SGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, WORK,
     $             RES )