sdrvsg2stg.f

Go to the documentation of this file.

*> \brief \b SDRVSG2STG
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SDRVSG2STG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
*                              NOUNIT, A, LDA, B, LDB, D, D2, Z, LDZ, AB,
*                              BB, AP, BP, WORK, NWORK, IWORK, LIWORK, 
*                              RESULT, INFO )
*
*       IMPLICIT NONE
*       .. Scalar Arguments ..
*       INTEGER            INFO, LDA, LDB, LDZ, LIWORK, NOUNIT, NSIZES,
*      $                   NTYPES, NWORK
*       REAL               THRESH
*       ..
*       .. Array Arguments ..
*       LOGICAL            DOTYPE( * )
*       INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
*       REAL               A( LDA, * ), AB( LDA, * ), AP( * ),
*      $                   B( LDB, * ), BB( LDB, * ), BP( * ), D( * ),
*      $                   RESULT( * ), WORK( * ), Z( LDZ, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>      SDRVSG2STG checks the real symmetric generalized eigenproblem
*>      drivers.
*>
*>              SSYGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem.
*>
*>              SSYGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem using a divide and conquer algorithm.
*>
*>              SSYGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem.
*>
*>              SSPGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage.
*>
*>              SSPGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage using a divide and
*>              conquer algorithm.
*>
*>              SSPGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite generalized
*>              eigenproblem in packed storage.
*>
*>              SSBGV computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem.
*>
*>              SSBGVD computes all eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem using a divide and conquer
*>              algorithm.
*>
*>              SSBGVX computes selected eigenvalues and, optionally,
*>              eigenvectors of a real symmetric-definite banded
*>              generalized eigenproblem.
*>
*>      When SDRVSG2STG is called, a number of matrix "sizes" ("n's") and a
*>      number of matrix "types" are specified.  For each size ("n")
*>      and each type of matrix, one matrix A of the given type will be
*>      generated; a random well-conditioned matrix B is also generated
*>      and the pair (A,B) is used to test the drivers.
*>
*>      For each pair (A,B), the following tests are performed:
*>
*>      (1) SSYGV with ITYPE = 1 and UPLO ='U':
*>
*>              | A Z - B Z D | / ( |A| |Z| n ulp )
*>              | D - D2 | / ( |D| ulp )   where D is computed by
*>                                         SSYGV and  D2 is computed by
*>                                         SSYGV_2STAGE. This test is
*>                                         only performed for SSYGV
*>
*>      (2) as (1) but calling SSPGV
*>      (3) as (1) but calling SSBGV
*>      (4) as (1) but with UPLO = 'L'
*>      (5) as (4) but calling SSPGV
*>      (6) as (4) but calling SSBGV
*>
*>      (7) SSYGV with ITYPE = 2 and UPLO ='U':
*>
*>              | A B Z - Z D | / ( |A| |Z| n ulp )
*>
*>      (8) as (7) but calling SSPGV
*>      (9) as (7) but with UPLO = 'L'
*>      (10) as (9) but calling SSPGV
*>
*>      (11) SSYGV with ITYPE = 3 and UPLO ='U':
*>
*>              | B A Z - Z D | / ( |A| |Z| n ulp )
*>
*>      (12) as (11) but calling SSPGV
*>      (13) as (11) but with UPLO = 'L'
*>      (14) as (13) but calling SSPGV
*>
*>      SSYGVD, SSPGVD and SSBGVD performed the same 14 tests.
*>
*>      SSYGVX, SSPGVX and SSBGVX performed the above 14 tests with
*>      the parameter RANGE = 'A', 'N' and 'I', respectively.
*>
*>      The "sizes" are specified by an array NN(1:NSIZES); the value
*>      of each element NN(j) specifies one size.
*>      The "types" are specified by a logical array DOTYPE( 1:NTYPES );
*>      if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
*>      This type is used for the matrix A which has half-bandwidth KA.
*>      B is generated as a well-conditioned positive definite matrix
*>      with half-bandwidth KB (<= KA).
*>      Currently, the list of possible types for A is:
*>
*>      (1)  The zero matrix.
*>      (2)  The identity matrix.
*>
*>      (3)  A diagonal matrix with evenly spaced entries
*>           1, ..., ULP  and random signs.
*>           (ULP = (first number larger than 1) - 1 )
*>      (4)  A diagonal matrix with geometrically spaced entries
*>           1, ..., ULP  and random signs.
*>      (5)  A diagonal matrix with "clustered" entries
*>           1, ULP, ..., ULP and random signs.
*>
*>      (6)  Same as (4), but multiplied by SQRT( overflow threshold )
*>      (7)  Same as (4), but multiplied by SQRT( underflow threshold )
*>
*>      (8)  A matrix of the form  U* D U, where U is orthogonal and
*>           D has evenly spaced entries 1, ..., ULP with random signs
*>           on the diagonal.
*>
*>      (9)  A matrix of the form  U* D U, where U is orthogonal and
*>           D has geometrically spaced entries 1, ..., ULP with random
*>           signs on the diagonal.
*>
*>      (10) A matrix of the form  U* D U, where U is orthogonal and
*>           D has "clustered" entries 1, ULP,..., ULP with random
*>           signs on the diagonal.
*>
*>      (11) Same as (8), but multiplied by SQRT( overflow threshold )
*>      (12) Same as (8), but multiplied by SQRT( underflow threshold )
*>
*>      (13) symmetric matrix with random entries chosen from (-1,1).
*>      (14) Same as (13), but multiplied by SQRT( overflow threshold )
*>      (15) Same as (13), but multiplied by SQRT( underflow threshold)
*>
*>      (16) Same as (8), but with KA = 1 and KB = 1
*>      (17) Same as (8), but with KA = 2 and KB = 1
*>      (18) Same as (8), but with KA = 2 and KB = 2
*>      (19) Same as (8), but with KA = 3 and KB = 1
*>      (20) Same as (8), but with KA = 3 and KB = 2
*>      (21) Same as (8), but with KA = 3 and KB = 3
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \verbatim
*>  NSIZES  INTEGER
*>          The number of sizes of matrices to use.  If it is zero,
*>          SDRVSG2STG does nothing.  It must be at least zero.
*>          Not modified.
*>
*>  NN      INTEGER array, dimension (NSIZES)
*>          An array containing the sizes to be used for the matrices.
*>          Zero values will be skipped.  The values must be at least
*>          zero.
*>          Not modified.
*>
*>  NTYPES  INTEGER
*>          The number of elements in DOTYPE.   If it is zero, SDRVSG2STG
*>          does nothing.  It must be at least zero.  If it is MAXTYP+1
*>          and NSIZES is 1, then an additional type, MAXTYP+1 is
*>          defined, which is to use whatever matrix is in A.  This
*>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
*>          DOTYPE(MAXTYP+1) is .TRUE. .
*>          Not modified.
*>
*>  DOTYPE  LOGICAL array, dimension (NTYPES)
*>          If DOTYPE(j) is .TRUE., then for each size in NN a
*>          matrix of that size and of type j will be generated.
*>          If NTYPES is smaller than the maximum number of types
*>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
*>          MAXTYP will not be generated.  If NTYPES is larger
*>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
*>          will be ignored.
*>          Not modified.
*>
*>  ISEED   INTEGER array, dimension (4)
*>          On entry ISEED specifies the seed of the random number
*>          generator. The array elements should be between 0 and 4095;
*>          if not they will be reduced mod 4096.  Also, ISEED(4) must
*>          be odd.  The random number generator uses a linear
*>          congruential sequence limited to small integers, and so
*>          should produce machine independent random numbers. The
*>          values of ISEED are changed on exit, and can be used in the
*>          next call to SDRVSG2STG to continue the same random number
*>          sequence.
*>          Modified.
*>
*>  THRESH  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. real)
*>          or the size of the matrix.  It must be at least zero.
*>          Not modified.
*>
*>  NOUNIT  INTEGER
*>          The FORTRAN unit number for printing out error messages
*>          (e.g., if a routine returns IINFO not equal to 0.)
*>          Not modified.
*>
*>  A       REAL             array, dimension (LDA , max(NN))
*>          Used to hold the matrix whose eigenvalues are to be
*>          computed.  On exit, A contains the last matrix actually
*>          used.
*>          Modified.
*>
*>  LDA     INTEGER
*>          The leading dimension of A and AB.  It must be at
*>          least 1 and at least max( NN ).
*>          Not modified.
*>
*>  B       REAL             array, dimension (LDB , max(NN))
*>          Used to hold the symmetric positive definite matrix for
*>          the generalized problem.
*>          On exit, B contains the last matrix actually
*>          used.
*>          Modified.
*>
*>  LDB     INTEGER
*>          The leading dimension of B and BB.  It must be at
*>          least 1 and at least max( NN ).
*>          Not modified.
*>
*>  D       REAL             array, dimension (max(NN))
*>          The eigenvalues of A. On exit, the eigenvalues in D
*>          correspond with the matrix in A.
*>          Modified.
*>
*>  Z       REAL             array, dimension (LDZ, max(NN))
*>          The matrix of eigenvectors.
*>          Modified.
*>
*>  LDZ     INTEGER
*>          The leading dimension of Z.  It must be at least 1 and
*>          at least max( NN ).
*>          Not modified.
*>
*>  AB      REAL             array, dimension (LDA, max(NN))
*>          Workspace.
*>          Modified.
*>
*>  BB      REAL             array, dimension (LDB, max(NN))
*>          Workspace.
*>          Modified.
*>
*>  AP      REAL             array, dimension (max(NN)**2)
*>          Workspace.
*>          Modified.
*>
*>  BP      REAL             array, dimension (max(NN)**2)
*>          Workspace.
*>          Modified.
*>
*>  WORK    REAL array, dimension (NWORK)
*>          Workspace.
*>          Modified.
*>
*>  NWORK   INTEGER
*>          The number of entries in WORK.  This must be at least
*>          1+5*N+2*N*lg(N)+3*N**2 where N = max( NN(j) ) and
*>          lg( N ) = smallest integer k such that 2**k >= N.
*>          Not modified.
*>
*>  IWORK   INTEGER array, dimension (LIWORK)
*>          Workspace.
*>          Modified.
*>
*>  LIWORK  INTEGER
*>          The number of entries in WORK.  This must be at least 6*N.
*>          Not modified.
*>
*>  RESULT  REAL array, dimension (70)
*>          The values computed by the 70 tests described above.
*>          Modified.
*>
*>  INFO    INTEGER
*>          If 0, then everything ran OK.
*>           -1: NSIZES < 0
*>           -2: Some NN(j) < 0
*>           -3: NTYPES < 0
*>           -5: THRESH < 0
*>           -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
*>          -16: LDZ < 1 or LDZ < NMAX.
*>          -21: NWORK too small.
*>          -23: LIWORK too small.
*>          If  SLATMR, SLATMS, SSYGV, SSPGV, SSBGV, SSYGVD, SSPGVD,
*>              SSBGVD, SSYGVX, SSPGVX or SSBGVX returns an error code,
*>              the absolute value of it is returned.
*>          Modified.
*>
*> ----------------------------------------------------------------------
*>
*>       Some Local Variables and Parameters:
*>       ---- ----- --------- --- ----------
*>       ZERO, ONE       Real 0 and 1.
*>       MAXTYP          The number of types defined.
*>       NTEST           The number of tests that have been run
*>                       on this matrix.
*>       NTESTT          The total number of tests for this call.
*>       NMAX            Largest value in NN.
*>       NMATS           The number of matrices generated so far.
*>       NERRS           The number of tests which have exceeded THRESH
*>                       so far (computed by SLAFTS).
*>       COND, IMODE     Values to be passed to the matrix generators.
*>       ANORM           Norm of A; passed to matrix generators.
*>
*>       OVFL, UNFL      Overflow and underflow thresholds.
*>       ULP, ULPINV     Finest relative precision and its inverse.
*>       RTOVFL, RTUNFL  Square roots of the previous 2 values.
*>               The following four arrays decode JTYPE:
*>       KTYPE(j)        The general type (1-10) for type "j".
*>       KMODE(j)        The MODE value to be passed to the matrix
*>                       generator for type "j".
*>       KMAGN(j)        The order of magnitude ( O(1),
*>                       O(overflow^(1/2) ), O(underflow^(1/2) )
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup real_eig
*
*  =====================================================================
      SUBROUTINE sdrvsg2stg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
     $                       NOUNIT, A, LDA, B, LDB, D, D2, Z, LDZ, AB,
     $                       BB, AP, BP, WORK, NWORK, IWORK, LIWORK, 
     $                       RESULT, INFO )
*
      IMPLICIT NONE
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            INFO, LDA, LDB, LDZ, LIWORK, NOUNIT, NSIZES,
     $                   NTYPES, NWORK
      REAL               THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            DOTYPE( * )
      INTEGER            ISEED( 4 ), IWORK( * ), NN( * )
      REAL               A( LDA, * ), AB( LDA, * ), AP( * ),
     $                   b( ldb, * ), bb( ldb, * ), bp( * ), d( * ),
     $                   d2( * ), result( * ), work( * ), z( ldz, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE, TEN
      PARAMETER          ( ZERO = 0.0e0, one = 1.0e0, ten = 10.0e0 )
      INTEGER            MAXTYP
      parameter( maxtyp = 21 )
*     ..
*     .. Local Scalars ..
      LOGICAL            BADNN
      CHARACTER          UPLO
      INTEGER            I, IBTYPE, IBUPLO, IINFO, IJ, IL, IMODE, ITEMP,
     $                   itype, iu, j, jcol, jsize, jtype, ka, ka9, kb,
     $                   kb9, m, mtypes, n, nerrs, nmats, nmax, ntest,
     $                   ntestt
      REAL               ABSTOL, ANINV, ANORM, COND, OVFL, RTOVFL,
     $                   RTUNFL, ULP, ULPINV, UNFL, VL, VU, TEMP1, TEMP2
*     ..
*     .. Local Arrays ..
      INTEGER            IDUMMA( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
     $                   KMAGN( MAXTYP ), KMODE( MAXTYP ),
     $                   KTYPE( MAXTYP )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      REAL               SLAMCH, SLARND
      EXTERNAL           LSAME, SLAMCH, SLARND
*     ..
*     .. External Subroutines ..
      EXTERNAL           slacpy, slafts, slaset, slasum, slatmr,
     $                   slatms, ssbgv, ssbgvd, ssbgvx, ssgt01, sspgv,
     $                   sspgvd, sspgvx, ssygv, ssygvd, ssygvx, xerbla,
     $                   ssygv_2stage
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, real, max, min, sqrt
*     ..
*     .. Data statements ..
      DATA               ktype / 1, 2, 5*4, 5*5, 3*8, 6*9 /
      DATA               kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
     $                   2, 3, 6*1 /
      DATA               kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
     $                   0, 0, 6*4 /
*     ..
*     .. Executable Statements ..
*
*     1)      Check for errors
*
      ntestt = 0
      info = 0
*
      badnn = .false.
      nmax = 0
      DO 10 j = 1, nsizes
         nmax = max( nmax, nn( j ) )
         IF( nn( j ).LT.0 )
     $      badnn = .true.
   10 CONTINUE
*
*     Check for errors
*
      IF( nsizes.LT.0 ) THEN
         info = -1
      ELSE IF( badnn ) THEN
         info = -2
      ELSE IF( ntypes.LT.0 ) THEN
         info = -3
      ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
         info = -9
      ELSE IF( ldz.LE.1 .OR. ldz.LT.nmax ) THEN
         info = -16
      ELSE IF( 2*max( nmax, 3 )**2.GT.nwork ) THEN
         info = -21
      ELSE IF( 2*max( nmax, 3 )**2.GT.liwork ) THEN
         info = -23
      END IF
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SDRVSG2STG', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
     $   RETURN
*
*     More Important constants
*
      unfl = slamch( 'Safe minimum' )
      ovfl = slamch( 'Overflow' )
      ulp = slamch( 'Epsilon' )*slamch( 'Base' )
      ulpinv = one / ulp
      rtunfl = sqrt( unfl )
      rtovfl = sqrt( ovfl )
*
      DO 20 i = 1, 4
         iseed2( i ) = iseed( i )
   20 CONTINUE
*
*     Loop over sizes, types
*
      nerrs = 0
      nmats = 0
*
      DO 650 jsize = 1, nsizes
         n = nn( jsize )
         aninv = one / real( max( 1, n ) )
*
         IF( nsizes.NE.1 ) THEN
            mtypes = min( maxtyp, ntypes )
         ELSE
            mtypes = min( maxtyp+1, ntypes )
         END IF
*
         ka9 = 0
         kb9 = 0
         DO 640 jtype = 1, mtypes
            IF( .NOT.dotype( jtype ) )
     $         GO TO 640
            nmats = nmats + 1
            ntest = 0
*
            DO 30 j = 1, 4
               ioldsd( j ) = iseed( j )
   30       CONTINUE
*
*           2)      Compute "A"
*
*                   Control parameters:
*
*               KMAGN  KMODE        KTYPE
*           =1  O(1)   clustered 1  zero
*           =2  large  clustered 2  identity
*           =3  small  exponential  (none)
*           =4         arithmetic   diagonal, w/ eigenvalues
*           =5         random log   hermitian, w/ eigenvalues
*           =6         random       (none)
*           =7                      random diagonal
*           =8                      random hermitian
*           =9                      banded, w/ eigenvalues
*
            IF( mtypes.GT.maxtyp )
     $         GO TO 90
*
            itype = ktype( jtype )
            imode = kmode( jtype )
*
*           Compute norm
*
            GO TO ( 40, 50, 60 )kmagn( jtype )
*
   40       CONTINUE
            anorm = one
            GO TO 70
*
   50       CONTINUE
            anorm = ( rtovfl*ulp )*aninv
            GO TO 70
*
   60       CONTINUE
            anorm = rtunfl*n*ulpinv
            GO TO 70
*
   70       CONTINUE
*
            iinfo = 0
            cond = ulpinv
*
*           Special Matrices -- Identity & Jordan block
*
            IF( itype.EQ.1 ) THEN
*
*              Zero
*
               ka = 0
               kb = 0
               CALL slaset( 'Full', lda, n, zero, zero, a, lda )
*
            ELSE IF( itype.EQ.2 ) THEN
*
*              Identity
*
               ka = 0
               kb = 0
               CALL slaset( 'Full', lda, n, zero, zero, a, lda )
               DO 80 jcol = 1, n
                  a( jcol, jcol ) = anorm
   80          CONTINUE
*
            ELSE IF( itype.EQ.4 ) THEN
*
*              Diagonal Matrix, [Eigen]values Specified
*
               ka = 0
               kb = 0
               CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
     $                      anorm, 0, 0, 'N', a, lda, work( n+1 ),
     $                      iinfo )
*
            ELSE IF( itype.EQ.5 ) THEN
*
*              symmetric, eigenvalues specified
*
               ka = max( 0, n-1 )
               kb = ka
               CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
     $                      anorm, n, n, 'N', a, lda, work( n+1 ),
     $                      iinfo )
*
            ELSE IF( itype.EQ.7 ) THEN
*
*              Diagonal, random eigenvalues
*
               ka = 0
               kb = 0
               CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
     $                      'T', 'N', work( n+1 ), 1, one,
     $                      work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
     $                      zero, anorm, 'NO', a, lda, iwork, iinfo )
*
            ELSE IF( itype.EQ.8 ) THEN
*
*              symmetric, random eigenvalues
*
               ka = max( 0, n-1 )
               kb = ka
               CALL slatmr( n, n, 'S', iseed, 'H', work, 6, one, one,
     $                      'T', 'N', work( n+1 ), 1, one,
     $                      work( 2*n+1 ), 1, one, 'N', idumma, n, n,
     $                      zero, anorm, 'NO', a, lda, iwork, iinfo )
*
            ELSE IF( itype.EQ.9 ) THEN
*
*              symmetric banded, eigenvalues specified
*
*              The following values are used for the half-bandwidths:
*
*                ka = 1   kb = 1
*                ka = 2   kb = 1
*                ka = 2   kb = 2
*                ka = 3   kb = 1
*                ka = 3   kb = 2
*                ka = 3   kb = 3
*
               kb9 = kb9 + 1
               IF( kb9.GT.ka9 ) THEN
                  ka9 = ka9 + 1
                  kb9 = 1
               END IF
               ka = max( 0, min( n-1, ka9 ) )
               kb = max( 0, min( n-1, kb9 ) )
               CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
     $                      anorm, ka, ka, 'N', a, lda, work( n+1 ),
     $                      iinfo )
*
            ELSE
*
               iinfo = 1
            END IF
*
            IF( iinfo.NE.0 ) THEN
               WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
     $            ioldsd
               info = abs( iinfo )
               RETURN
            END IF
*
   90       CONTINUE
*
            abstol = unfl + unfl
            IF( n.LE.1 ) THEN
               il = 1
               iu = n
            ELSE
               il = 1 + int( ( n-1 )*slarnd( 1, iseed2 ) )
               iu = 1 + int( ( n-1 )*slarnd( 1, iseed2 ) )
               IF( il.GT.iu ) THEN
                  itemp = il
                  il = iu
                  iu = itemp
               END IF
            END IF
*
*           3) Call SSYGV, SSPGV, SSBGV, SSYGVD, SSPGVD, SSBGVD,
*              SSYGVX, SSPGVX, and SSBGVX, do tests.
*
*           loop over the three generalized problems
*                 IBTYPE = 1: A*x = (lambda)*B*x
*                 IBTYPE = 2: A*B*x = (lambda)*x
*                 IBTYPE = 3: B*A*x = (lambda)*x
*
            DO 630 ibtype = 1, 3
*
*              loop over the setting UPLO
*
               DO 620 ibuplo = 1, 2
                  IF( ibuplo.EQ.1 )
     $               uplo = 'U'
                  IF( ibuplo.EQ.2 )
     $               uplo = 'L'
*
*                 Generate random well-conditioned positive definite
*                 matrix B, of bandwidth not greater than that of A.
*
                  CALL slatms( n, n, 'U', iseed, 'P', work, 5, ten, one,
     $                         kb, kb, uplo, b, ldb, work( n+1 ),
     $                         iinfo )
*
*                 Test SSYGV
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, z, ldz )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
                  CALL ssygv( ibtype, 'V', uplo, n, z, ldz, bb, ldb, d,
     $                        work, nwork, iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSYGV(V,' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
*                 Test SSYGV_2STAGE
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, z, ldz )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
                  CALL ssygv_2stage( ibtype, 'N', uplo, n, z, ldz,
     $                               bb, ldb, d2, work, nwork, iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )
     $                  'SSYGV_2STAGE(V,' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
C                  CALL SSGT01( IBTYPE, UPLO, N, N, A, LDA, B, LDB, Z,
C     $                         LDZ, D, WORK, RESULT( NTEST ) )
*           
*           
*                 Do Tests | D1 - D2 | / ( |D1| ulp )
*                 D1 computed using the standard 1-stage reduction as reference
*                 D2 computed using the 2-stage reduction
*           
                  temp1 = zero
                  temp2 = zero
                  DO 151 j = 1, n
                     temp1 = max( temp1, abs( d( j ) ), 
     $                                   abs( d2( j ) ) )
                     temp2 = max( temp2, abs( d( j )-d2( j ) ) )
  151             CONTINUE
*           
                  result( ntest ) = temp2 / 
     $                              max( unfl, ulp*max( temp1, temp2 ) )
*
*                 Test SSYGVD
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, z, ldz )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
                  CALL ssygvd( ibtype, 'V', uplo, n, z, ldz, bb, ldb, d,
     $                         work, nwork, iwork, liwork, iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSYGVD(V,' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
*                 Test SSYGVX
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, ab, lda )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
                  CALL ssygvx( ibtype, 'V', 'A', uplo, n, ab, lda, bb,
     $                         ldb, vl, vu, il, iu, abstol, m, d, z,
     $                         ldz, work, nwork, iwork( n+1 ), iwork,
     $                         iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSYGVX(V,A' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, ab, lda )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
*                 since we do not know the exact eigenvalues of this
*                 eigenpair, we just set VL and VU as constants.
*                 It is quite possible that there are no eigenvalues
*                 in this interval.
*
                  vl = zero
                  vu = anorm
                  CALL ssygvx( ibtype, 'V', 'V', uplo, n, ab, lda, bb,
     $                         ldb, vl, vu, il, iu, abstol, m, d, z,
     $                         ldz, work, nwork, iwork( n+1 ), iwork,
     $                         iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSYGVX(V,V,' //
     $                  uplo // ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
                  ntest = ntest + 1
*
                  CALL slacpy( ' ', n, n, a, lda, ab, lda )
                  CALL slacpy( uplo, n, n, b, ldb, bb, ldb )
*
                  CALL ssygvx( ibtype, 'V', 'I', uplo, n, ab, lda, bb,
     $                         ldb, vl, vu, il, iu, abstol, m, d, z,
     $                         ldz, work, nwork, iwork( n+1 ), iwork,
     $                         iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSYGVX(V,I,' //
     $                  uplo // ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 100
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
  100             CONTINUE
*
*                 Test SSPGV
*
                  ntest = ntest + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( lsame( uplo, 'U' ) ) THEN
                     ij = 1
                     DO 120 j = 1, n
                        DO 110 i = 1, j
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  110                   CONTINUE
  120                CONTINUE
                  ELSE
                     ij = 1
                     DO 140 j = 1, n
                        DO 130 i = j, n
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  130                   CONTINUE
  140                CONTINUE
                  END IF
*
                  CALL sspgv( ibtype, 'V', uplo, n, ap, bp, d, z, ldz,
     $                        work, iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSPGV(V,' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
*                 Test SSPGVD
*
                  ntest = ntest + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( lsame( uplo, 'U' ) ) THEN
                     ij = 1
                     DO 160 j = 1, n
                        DO 150 i = 1, j
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  150                   CONTINUE
  160                CONTINUE
                  ELSE
                     ij = 1
                     DO 180 j = 1, n
                        DO 170 i = j, n
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  170                   CONTINUE
  180                CONTINUE
                  END IF
*
                  CALL sspgvd( ibtype, 'V', uplo, n, ap, bp, d, z, ldz,
     $                         work, nwork, iwork, liwork, iinfo )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSPGVD(V,' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
*                 Test SSPGVX
*
                  ntest = ntest + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( lsame( uplo, 'U' ) ) THEN
                     ij = 1
                     DO 200 j = 1, n
                        DO 190 i = 1, j
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  190                   CONTINUE
  200                CONTINUE
                  ELSE
                     ij = 1
                     DO 220 j = 1, n
                        DO 210 i = j, n
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  210                   CONTINUE
  220                CONTINUE
                  END IF
*
                  CALL sspgvx( ibtype, 'V', 'A', uplo, n, ap, bp, vl,
     $                         vu, il, iu, abstol, m, d, z, ldz, work,
     $                         iwork( n+1 ), iwork, info )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSPGVX(V,A' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
                  ntest = ntest + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( lsame( uplo, 'U' ) ) THEN
                     ij = 1
                     DO 240 j = 1, n
                        DO 230 i = 1, j
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  230                   CONTINUE
  240                CONTINUE
                  ELSE
                     ij = 1
                     DO 260 j = 1, n
                        DO 250 i = j, n
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  250                   CONTINUE
  260                CONTINUE
                  END IF
*
                  vl = zero
                  vu = anorm
                  CALL sspgvx( ibtype, 'V', 'V', uplo, n, ap, bp, vl,
     $                         vu, il, iu, abstol, m, d, z, ldz, work,
     $                         iwork( n+1 ), iwork, info )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSPGVX(V,V' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
                  ntest = ntest + 1
*
*                 Copy the matrices into packed storage.
*
                  IF( lsame( uplo, 'U' ) ) THEN
                     ij = 1
                     DO 280 j = 1, n
                        DO 270 i = 1, j
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  270                   CONTINUE
  280                CONTINUE
                  ELSE
                     ij = 1
                     DO 300 j = 1, n
                        DO 290 i = j, n
                           ap( ij ) = a( i, j )
                           bp( ij ) = b( i, j )
                           ij = ij + 1
  290                   CONTINUE
  300                CONTINUE
                  END IF
*
                  CALL sspgvx( ibtype, 'V', 'I', uplo, n, ap, bp, vl,
     $                         vu, il, iu, abstol, m, d, z, ldz, work,
     $                         iwork( n+1 ), iwork, info )
                  IF( iinfo.NE.0 ) THEN
                     WRITE( nounit, fmt = 9999 )'SSPGVX(V,I' // uplo //
     $                  ')', iinfo, n, jtype, ioldsd
                     info = abs( iinfo )
                     IF( iinfo.LT.0 ) THEN
                        RETURN
                     ELSE
                        result( ntest ) = ulpinv
                        GO TO 310
                     END IF
                  END IF
*
*                 Do Test
*
                  CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                         ldz, d, work, result( ntest ) )
*
  310             CONTINUE
*
                  IF( ibtype.EQ.1 ) THEN
*
*                    TEST SSBGV
*
                     ntest = ntest + 1
*
*                    Copy the matrices into band storage.
*
                     IF( lsame( uplo, 'U' ) ) THEN
                        DO 340 j = 1, n
                           DO 320 i = max( 1, j-ka ), j
                              ab( ka+1+i-j, j ) = a( i, j )
  320                      CONTINUE
                           DO 330 i = max( 1, j-kb ), j
                              bb( kb+1+i-j, j ) = b( i, j )
  330                      CONTINUE
  340                   CONTINUE
                     ELSE
                        DO 370 j = 1, n
                           DO 350 i = j, min( n, j+ka )
                              ab( 1+i-j, j ) = a( i, j )
  350                      CONTINUE
                           DO 360 i = j, min( n, j+kb )
                              bb( 1+i-j, j ) = b( i, j )
  360                      CONTINUE
  370                   CONTINUE
                     END IF
*
                     CALL ssbgv( 'V', uplo, n, ka, kb, ab, lda, bb, ldb,
     $                           d, z, ldz, work, iinfo )
                     IF( iinfo.NE.0 ) THEN
                        WRITE( nounit, fmt = 9999 )'SSBGV(V,' //
     $                     uplo // ')', iinfo, n, jtype, ioldsd
                        info = abs( iinfo )
                        IF( iinfo.LT.0 ) THEN
                           RETURN
                        ELSE
                           result( ntest ) = ulpinv
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                            ldz, d, work, result( ntest ) )
*
*                    TEST SSBGVD
*
                     ntest = ntest + 1
*
*                    Copy the matrices into band storage.
*
                     IF( lsame( uplo, 'U' ) ) THEN
                        DO 400 j = 1, n
                           DO 380 i = max( 1, j-ka ), j
                              ab( ka+1+i-j, j ) = a( i, j )
  380                      CONTINUE
                           DO 390 i = max( 1, j-kb ), j
                              bb( kb+1+i-j, j ) = b( i, j )
  390                      CONTINUE
  400                   CONTINUE
                     ELSE
                        DO 430 j = 1, n
                           DO 410 i = j, min( n, j+ka )
                              ab( 1+i-j, j ) = a( i, j )
  410                      CONTINUE
                           DO 420 i = j, min( n, j+kb )
                              bb( 1+i-j, j ) = b( i, j )
  420                      CONTINUE
  430                   CONTINUE
                     END IF
*
                     CALL ssbgvd( 'V', uplo, n, ka, kb, ab, lda, bb,
     $                            ldb, d, z, ldz, work, nwork, iwork,
     $                            liwork, iinfo )
                     IF( iinfo.NE.0 ) THEN
                        WRITE( nounit, fmt = 9999 )'SSBGVD(V,' //
     $                     uplo // ')', iinfo, n, jtype, ioldsd
                        info = abs( iinfo )
                        IF( iinfo.LT.0 ) THEN
                           RETURN
                        ELSE
                           result( ntest ) = ulpinv
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL ssgt01( ibtype, uplo, n, n, a, lda, b, ldb, z,
     $                            ldz, d, work, result( ntest ) )
*
*                    Test SSBGVX
*
                     ntest = ntest + 1
*
*                    Copy the matrices into band storage.
*
                     IF( lsame( uplo, 'U' ) ) THEN
                        DO 460 j = 1, n
                           DO 440 i = max( 1, j-ka ), j
                              ab( ka+1+i-j, j ) = a( i, j )
  440                      CONTINUE
                           DO 450 i = max( 1, j-kb ), j
                              bb( kb+1+i-j, j ) = b( i, j )
  450                      CONTINUE
  460                   CONTINUE
                     ELSE
                        DO 490 j = 1, n
                           DO 470 i = j, min( n, j+ka )
                              ab( 1+i-j, j ) = a( i, j )
  470                      CONTINUE
                           DO 480 i = j, min( n, j+kb )
                              bb( 1+i-j, j ) = b( i, j )
  480                      CONTINUE
  490                   CONTINUE
                     END IF
*
                     CALL ssbgvx( 'V', 'A', uplo, n, ka, kb, ab, lda,
     $                            bb, ldb, bp, max( 1, n ), vl, vu, il,
     $                            iu, abstol, m, d, z, ldz, work,
     $                            iwork( n+1 ), iwork, iinfo )
                     IF( iinfo.NE.0 ) THEN
                        WRITE( nounit, fmt = 9999 )'SSBGVX(V,A' //
     $                     uplo // ')', iinfo, n, jtype, ioldsd
                        info = abs( iinfo )
                        IF( iinfo.LT.0 ) THEN
                           RETURN
                        ELSE
                           result( ntest ) = ulpinv
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                            ldz, d, work, result( ntest ) )
*
*
                     ntest = ntest + 1
*
*                    Copy the matrices into band storage.
*
                     IF( lsame( uplo, 'U' ) ) THEN
                        DO 520 j = 1, n
                           DO 500 i = max( 1, j-ka ), j
                              ab( ka+1+i-j, j ) = a( i, j )
  500                      CONTINUE
                           DO 510 i = max( 1, j-kb ), j
                              bb( kb+1+i-j, j ) = b( i, j )
  510                      CONTINUE
  520                   CONTINUE
                     ELSE
                        DO 550 j = 1, n
                           DO 530 i = j, min( n, j+ka )
                              ab( 1+i-j, j ) = a( i, j )
  530                      CONTINUE
                           DO 540 i = j, min( n, j+kb )
                              bb( 1+i-j, j ) = b( i, j )
  540                      CONTINUE
  550                   CONTINUE
                     END IF
*
                     vl = zero
                     vu = anorm
                     CALL ssbgvx( 'V', 'V', uplo, n, ka, kb, ab, lda,
     $                            bb, ldb, bp, max( 1, n ), vl, vu, il,
     $                            iu, abstol, m, d, z, ldz, work,
     $                            iwork( n+1 ), iwork, iinfo )
                     IF( iinfo.NE.0 ) THEN
                        WRITE( nounit, fmt = 9999 )'SSBGVX(V,V' //
     $                     uplo // ')', iinfo, n, jtype, ioldsd
                        info = abs( iinfo )
                        IF( iinfo.LT.0 ) THEN
                           RETURN
                        ELSE
                           result( ntest ) = ulpinv
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                            ldz, d, work, result( ntest ) )
*
                     ntest = ntest + 1
*
*                    Copy the matrices into band storage.
*
                     IF( lsame( uplo, 'U' ) ) THEN
                        DO 580 j = 1, n
                           DO 560 i = max( 1, j-ka ), j
                              ab( ka+1+i-j, j ) = a( i, j )
  560                      CONTINUE
                           DO 570 i = max( 1, j-kb ), j
                              bb( kb+1+i-j, j ) = b( i, j )
  570                      CONTINUE
  580                   CONTINUE
                     ELSE
                        DO 610 j = 1, n
                           DO 590 i = j, min( n, j+ka )
                              ab( 1+i-j, j ) = a( i, j )
  590                      CONTINUE
                           DO 600 i = j, min( n, j+kb )
                              bb( 1+i-j, j ) = b( i, j )
  600                      CONTINUE
  610                   CONTINUE
                     END IF
*
                     CALL ssbgvx( 'V', 'I', uplo, n, ka, kb, ab, lda,
     $                            bb, ldb, bp, max( 1, n ), vl, vu, il,
     $                            iu, abstol, m, d, z, ldz, work,
     $                            iwork( n+1 ), iwork, iinfo )
                     IF( iinfo.NE.0 ) THEN
                        WRITE( nounit, fmt = 9999 )'SSBGVX(V,I' //
     $                     uplo // ')', iinfo, n, jtype, ioldsd
                        info = abs( iinfo )
                        IF( iinfo.LT.0 ) THEN
                           RETURN
                        ELSE
                           result( ntest ) = ulpinv
                           GO TO 620
                        END IF
                     END IF
*
*                    Do Test
*
                     CALL ssgt01( ibtype, uplo, n, m, a, lda, b, ldb, z,
     $                            ldz, d, work, result( ntest ) )
*
                  END IF
*
  620          CONTINUE
  630       CONTINUE
*
*           End of Loop -- Check for RESULT(j) > THRESH
*
            ntestt = ntestt + ntest
            CALL slafts( 'SSG', n, n, jtype, ntest, result, ioldsd,
     $                   thresh, nounit, nerrs )
  640    CONTINUE
  650 CONTINUE
*
*     Summary
*
      CALL slasum( 'SSG', nounit, nerrs, ntestt )
*
      RETURN
*
*     End of SDRVSG2STG
*
 9999 FORMAT( ' SDRVSG2STG: ', a, ' returned INFO=', i6, '.', / 9x,
     $    'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
      END