LAPACK: TESTING/EIG/sdrges.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SDRGES( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00002      $                   NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR,
00003      $                   ALPHAI, BETA, WORK, LWORK, RESULT, BWORK,
00004      $                   INFO )
00005 *
00006 *  -- LAPACK test routine (version 3.1) --
00007 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       INTEGER            INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
00012       REAL               THRESH
00013 *     ..
00014 *     .. Array Arguments ..
00015       LOGICAL            BWORK( * ), DOTYPE( * )
00016       INTEGER            ISEED( 4 ), NN( * )
00017       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00018      $                   B( LDA, * ), BETA( * ), Q( LDQ, * ),
00019      $                   RESULT( 13 ), S( LDA, * ), T( LDA, * ),
00020      $                   WORK( * ), Z( LDQ, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  SDRGES checks the nonsymmetric generalized eigenvalue (Schur form)
00027 *  problem driver SGGES.
00028 *
00029 *  SGGES factors A and B as Q S Z'  and Q T Z' , where ' means
00030 *  transpose, T is upper triangular, S is in generalized Schur form
00031 *  (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
00032 *  the 2x2 blocks corresponding to complex conjugate pairs of
00033 *  generalized eigenvalues), and Q and Z are orthogonal. It also
00034 *  computes the generalized eigenvalues (alpha(j),beta(j)), j=1,...,n,
00035 *  Thus, w(j) = alpha(j)/beta(j) is a root of the characteristic
00036 *  equation
00037 *                  det( A - w(j) B ) = 0
00038 *  Optionally it also reorder the eigenvalues so that a selected
00039 *  cluster of eigenvalues appears in the leading diagonal block of the
00040 *  Schur forms.
00041 *
00042 *  When SDRGES is called, a number of matrix "sizes" ("N's") and a
00043 *  number of matrix "TYPES" are specified.  For each size ("N")
00044 *  and each TYPE of matrix, a pair of matrices (A, B) will be generated
00045 *  and used for testing. For each matrix pair, the following 13 tests
00046 *  will be performed and compared with the threshhold THRESH except
00047 *  the tests (5), (11) and (13).
00048 *
00049 *
00050 *  (1)   | A - Q S Z' | / ( |A| n ulp ) (no sorting of eigenvalues)
00051 *
00052 *
00053 *  (2)   | B - Q T Z' | / ( |B| n ulp ) (no sorting of eigenvalues)
00054 *
00055 *
00056 *  (3)   | I - QQ' | / ( n ulp ) (no sorting of eigenvalues)
00057 *
00058 *
00059 *  (4)   | I - ZZ' | / ( n ulp ) (no sorting of eigenvalues)
00060 *
00061 *  (5)   if A is in Schur form (i.e. quasi-triangular form)
00062 *        (no sorting of eigenvalues)
00063 *
00064 *  (6)   if eigenvalues = diagonal blocks of the Schur form (S, T),
00065 *        i.e., test the maximum over j of D(j)  where:
00066 *
00067 *        if alpha(j) is real:
00068 *                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
00069 *            D(j) = ------------------------ + -----------------------
00070 *                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
00071 *
00072 *        if alpha(j) is complex:
00073 *                                  | det( s S - w T ) |
00074 *            D(j) = ---------------------------------------------------
00075 *                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
00076 *
00077 *        and S and T are here the 2 x 2 diagonal blocks of S and T
00078 *        corresponding to the j-th and j+1-th eigenvalues.
00079 *        (no sorting of eigenvalues)
00080 *
00081 *  (7)   | (A,B) - Q (S,T) Z' | / ( | (A,B) | n ulp )
00082 *             (with sorting of eigenvalues).
00083 *
00084 *  (8)   | I - QQ' | / ( n ulp ) (with sorting of eigenvalues).
00085 *
00086 *  (9)   | I - ZZ' | / ( n ulp ) (with sorting of eigenvalues).
00087 *
00088 *  (10)  if A is in Schur form (i.e. quasi-triangular form)
00089 *        (with sorting of eigenvalues).
00090 *
00091 *  (11)  if eigenvalues = diagonal blocks of the Schur form (S, T),
00092 *        i.e. test the maximum over j of D(j)  where:
00093 *
00094 *        if alpha(j) is real:
00095 *                      |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
00096 *            D(j) = ------------------------ + -----------------------
00097 *                   max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
00098 *
00099 *        if alpha(j) is complex:
00100 *                                  | det( s S - w T ) |
00101 *            D(j) = ---------------------------------------------------
00102 *                   ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
00103 *
00104 *        and S and T are here the 2 x 2 diagonal blocks of S and T
00105 *        corresponding to the j-th and j+1-th eigenvalues.
00106 *        (with sorting of eigenvalues).
00107 *
00108 *  (12)  if sorting worked and SDIM is the number of eigenvalues
00109 *        which were SELECTed.
00110 *
00111 *  Test Matrices
00112 *  =============
00113 *
00114 *  The sizes of the test matrices are specified by an array
00115 *  NN(1:NSIZES); the value of each element NN(j) specifies one size.
00116 *  The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
00117 *  DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00118 *  Currently, the list of possible types is:
00119 *
00120 *  (1)  ( 0, 0 )         (a pair of zero matrices)
00121 *
00122 *  (2)  ( I, 0 )         (an identity and a zero matrix)
00123 *
00124 *  (3)  ( 0, I )         (an identity and a zero matrix)
00125 *
00126 *  (4)  ( I, I )         (a pair of identity matrices)
00127 *
00128 *          t   t
00129 *  (5)  ( J , J  )       (a pair of transposed Jordan blocks)
00130 *
00131 *                                      t                ( I   0  )
00132 *  (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t )
00133 *                                   ( 0   I  )          ( 0   J  )
00134 *                        and I is a k x k identity and J a (k+1)x(k+1)
00135 *                        Jordan block; k=(N-1)/2
00136 *
00137 *  (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal
00138 *                        matrix with those diagonal entries.)
00139 *  (8)  ( I, D )
00140 *
00141 *  (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big
00142 *
00143 *  (10) ( small*D, big*I )
00144 *
00145 *  (11) ( big*I, small*D )
00146 *
00147 *  (12) ( small*I, big*D )
00148 *
00149 *  (13) ( big*D, big*I )
00150 *
00151 *  (14) ( small*D, small*I )
00152 *
00153 *  (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
00154 *                         D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
00155 *            t   t
00156 *  (16) Q ( J , J ) Z     where Q and Z are random orthogonal matrices.
00157 *
00158 *  (17) Q ( T1, T2 ) Z    where T1 and T2 are upper triangular matrices
00159 *                         with random O(1) entries above the diagonal
00160 *                         and diagonal entries diag(T1) =
00161 *                         ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
00162 *                         ( 0, N-3, N-4,..., 1, 0, 0 )
00163 *
00164 *  (18) Q ( T1, T2 ) Z    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
00165 *                         diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
00166 *                         s = machine precision.
00167 *
00168 *  (19) Q ( T1, T2 ) Z    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
00169 *                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
00170 *
00171 *                                                         N-5
00172 *  (20) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 )
00173 *                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
00174 *
00175 *  (21) Q ( T1, T2 ) Z    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
00176 *                         diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
00177 *                         where r1,..., r(N-4) are random.
00178 *
00179 *  (22) Q ( big*T1, small*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00180 *                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00181 *
00182 *  (23) Q ( small*T1, big*T2 ) Z    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00183 *                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00184 *
00185 *  (24) Q ( small*T1, small*T2 ) Z  diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00186 *                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00187 *
00188 *  (25) Q ( big*T1, big*T2 ) Z      diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00189 *                                   diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00190 *
00191 *  (26) Q ( T1, T2 ) Z     where T1 and T2 are random upper-triangular
00192 *                          matrices.
00193 *
00194 *
00195 *  Arguments
00196 *  =========
00197 *
00198 *  NSIZES  (input) INTEGER
00199 *          The number of sizes of matrices to use.  If it is zero,
00200 *          SDRGES does nothing.  NSIZES >= 0.
00201 *
00202 *  NN      (input) INTEGER array, dimension (NSIZES)
00203 *          An array containing the sizes to be used for the matrices.
00204 *          Zero values will be skipped.  NN >= 0.
00205 *
00206 *  NTYPES  (input) INTEGER
00207 *          The number of elements in DOTYPE.   If it is zero, SDRGES
00208 *          does nothing.  It must be at least zero.  If it is MAXTYP+1
00209 *          and NSIZES is 1, then an additional type, MAXTYP+1 is
00210 *          defined, which is to use whatever matrix is in A on input.
00211 *          This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00212 *          DOTYPE(MAXTYP+1) is .TRUE. .
00213 *
00214 *  DOTYPE  (input) LOGICAL array, dimension (NTYPES)
00215 *          If DOTYPE(j) is .TRUE., then for each size in NN a
00216 *          matrix of that size and of type j will be generated.
00217 *          If NTYPES is smaller than the maximum number of types
00218 *          defined (PARAMETER MAXTYP), then types NTYPES+1 through
00219 *          MAXTYP will not be generated. If NTYPES is larger
00220 *          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00221 *          will be ignored.
00222 *
00223 *  ISEED   (input/output) INTEGER array, dimension (4)
00224 *          On entry ISEED specifies the seed of the random number
00225 *          generator. The array elements should be between 0 and 4095;
00226 *          if not they will be reduced mod 4096. Also, ISEED(4) must
00227 *          be odd.  The random number generator uses a linear
00228 *          congruential sequence limited to small integers, and so
00229 *          should produce machine independent random numbers. The
00230 *          values of ISEED are changed on exit, and can be used in the
00231 *          next call to SDRGES to continue the same random number
00232 *          sequence.
00233 *
00234 *  THRESH  (input) REAL
00235 *          A test will count as "failed" if the "error", computed as
00236 *          described above, exceeds THRESH.  Note that the error is
00237 *          scaled to be O(1), so THRESH should be a reasonably small
00238 *          multiple of 1, e.g., 10 or 100.  In particular, it should
00239 *          not depend on the precision (single vs. double) or the size
00240 *          of the matrix.  THRESH >= 0.
00241 *
00242 *  NOUNIT  (input) INTEGER
00243 *          The FORTRAN unit number for printing out error messages
00244 *          (e.g., if a routine returns IINFO not equal to 0.)
00245 *
00246 *  A       (input/workspace) REAL array,
00247 *                                       dimension(LDA, max(NN))
00248 *          Used to hold the original A matrix.  Used as input only
00249 *          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
00250 *          DOTYPE(MAXTYP+1)=.TRUE.
00251 *
00252 *  LDA     (input) INTEGER
00253 *          The leading dimension of A, B, S, and T.
00254 *          It must be at least 1 and at least max( NN ).
00255 *
00256 *  B       (input/workspace) REAL array,
00257 *                                       dimension(LDA, max(NN))
00258 *          Used to hold the original B matrix.  Used as input only
00259 *          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
00260 *          DOTYPE(MAXTYP+1)=.TRUE.
00261 *
00262 *  S       (workspace) REAL array, dimension (LDA, max(NN))
00263 *          The Schur form matrix computed from A by SGGES.  On exit, S
00264 *          contains the Schur form matrix corresponding to the matrix
00265 *          in A.
00266 *
00267 *  T       (workspace) REAL array, dimension (LDA, max(NN))
00268 *          The upper triangular matrix computed from B by SGGES.
00269 *
00270 *  Q       (workspace) REAL array, dimension (LDQ, max(NN))
00271 *          The (left) orthogonal matrix computed by SGGES.
00272 *
00273 *  LDQ     (input) INTEGER
00274 *          The leading dimension of Q and Z. It must
00275 *          be at least 1 and at least max( NN ).
00276 *
00277 *  Z       (workspace) REAL array, dimension( LDQ, max(NN) )
00278 *          The (right) orthogonal matrix computed by SGGES.
00279 *
00280 *  ALPHAR  (workspace) REAL array, dimension (max(NN))
00281 *  ALPHAI  (workspace) REAL array, dimension (max(NN))
00282 *  BETA    (workspace) REAL array, dimension (max(NN))
00283 *          The generalized eigenvalues of (A,B) computed by SGGES.
00284 *          ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
00285 *          generalized eigenvalue of A and B.
00286 *
00287 *  WORK    (workspace) REAL array, dimension (LWORK)
00288 *
00289 *  LWORK   (input) INTEGER
00290 *          The dimension of the array WORK.
00291 *          LWORK >= MAX( 10*(N+1), 3*N*N ), where N is the largest
00292 *          matrix dimension.
00293 *
00294 *  RESULT  (output) REAL array, dimension (15)
00295 *          The values computed by the tests described above.
00296 *          The values are currently limited to 1/ulp, to avoid overflow.
00297 *
00298 *  BWORK   (workspace) LOGICAL array, dimension (N)
00299 *
00300 *  INFO    (output) INTEGER
00301 *          = 0:  successful exit
00302 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00303 *          > 0:  A routine returned an error code.  INFO is the
00304 *                absolute value of the INFO value returned.
00305 *
00306 *  =====================================================================
00307 *
00308 *     .. Parameters ..
00309       REAL               ZERO, ONE
00310       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00311       INTEGER            MAXTYP
00312       PARAMETER          ( MAXTYP = 26 )
00313 *     ..
00314 *     .. Local Scalars ..
00315       LOGICAL            BADNN, ILABAD
00316       CHARACTER          SORT
00317       INTEGER            I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
00318      $                   JSIZE, JTYPE, KNTEIG, MAXWRK, MINWRK, MTYPES,
00319      $                   N, N1, NB, NERRS, NMATS, NMAX, NTEST, NTESTT,
00320      $                   RSUB, SDIM
00321       REAL               SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
00322 *     ..
00323 *     .. Local Arrays ..
00324       INTEGER            IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
00325      $                   IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
00326      $                   KATYPE( MAXTYP ), KAZERO( MAXTYP ),
00327      $                   KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
00328      $                   KBZERO( MAXTYP ), KCLASS( MAXTYP ),
00329      $                   KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
00330       REAL               RMAGN( 0: 3 )
00331 *     ..
00332 *     .. External Functions ..
00333       LOGICAL            SLCTES
00334       INTEGER            ILAENV
00335       REAL               SLAMCH, SLARND
00336       EXTERNAL           SLCTES, ILAENV, SLAMCH, SLARND
00337 *     ..
00338 *     .. External Subroutines ..
00339       EXTERNAL           ALASVM, SGET51, SGET53, SGET54, SGGES, SLABAD,
00340      $                   SLACPY, SLARFG, SLASET, SLATM4, SORM2R, XERBLA
00341 *     ..
00342 *     .. Intrinsic Functions ..
00343       INTRINSIC          ABS, MAX, MIN, REAL, SIGN
00344 *     ..
00345 *     .. Data statements ..
00346       DATA               KCLASS / 15*1, 10*2, 1*3 /
00347       DATA               KZ1 / 0, 1, 2, 1, 3, 3 /
00348       DATA               KZ2 / 0, 0, 1, 2, 1, 1 /
00349       DATA               KADD / 0, 0, 0, 0, 3, 2 /
00350       DATA               KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
00351      $                   4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
00352       DATA               KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
00353      $                   1, 1, -4, 2, -4, 8*8, 0 /
00354       DATA               KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
00355      $                   4*5, 4*3, 1 /
00356       DATA               KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
00357      $                   4*6, 4*4, 1 /
00358       DATA               KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
00359      $                   2, 1 /
00360       DATA               KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
00361      $                   2, 1 /
00362       DATA               KTRIAN / 16*0, 10*1 /
00363       DATA               IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
00364      $                   5*2, 0 /
00365       DATA               IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
00366 *     ..
00367 *     .. Executable Statements ..
00368 *
00369 *     Check for errors
00370 *
00371       INFO = 0
00372 *
00373       BADNN = .FALSE.
00374       NMAX = 1
00375       DO 10 J = 1, NSIZES
00376          NMAX = MAX( NMAX, NN( J ) )
00377          IF( NN( J ).LT.0 )
00378      $      BADNN = .TRUE.
00379    10 CONTINUE
00380 *
00381       IF( NSIZES.LT.0 ) THEN
00382          INFO = -1
00383       ELSE IF( BADNN ) THEN
00384          INFO = -2
00385       ELSE IF( NTYPES.LT.0 ) THEN
00386          INFO = -3
00387       ELSE IF( THRESH.LT.ZERO ) THEN
00388          INFO = -6
00389       ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
00390          INFO = -9
00391       ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
00392          INFO = -14
00393       END IF
00394 *
00395 *     Compute workspace
00396 *      (Note: Comments in the code beginning "Workspace:" describe the
00397 *       minimal amount of workspace needed at that point in the code,
00398 *       as well as the preferred amount for good performance.
00399 *       NB refers to the optimal block size for the immediately
00400 *       following subroutine, as returned by ILAENV.
00401 *
00402       MINWRK = 1
00403       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00404          MINWRK = MAX( 10*( NMAX+1 ), 3*NMAX*NMAX )
00405          NB = MAX( 1, ILAENV( 1, 'SGEQRF', ' ', NMAX, NMAX, -1, -1 ),
00406      $        ILAENV( 1, 'SORMQR', 'LT', NMAX, NMAX, NMAX, -1 ),
00407      $        ILAENV( 1, 'SORGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
00408          MAXWRK = MAX( 10*( NMAX+1 ), 2*NMAX+NMAX*NB, 3*NMAX*NMAX )
00409          WORK( 1 ) = MAXWRK
00410       END IF
00411 *
00412       IF( LWORK.LT.MINWRK )
00413      $   INFO = -20
00414 *
00415       IF( INFO.NE.0 ) THEN
00416          CALL XERBLA( 'SDRGES', -INFO )
00417          RETURN
00418       END IF
00419 *
00420 *     Quick return if possible
00421 *
00422       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00423      $   RETURN
00424 *
00425       SAFMIN = SLAMCH( 'Safe minimum' )
00426       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00427       SAFMIN = SAFMIN / ULP
00428       SAFMAX = ONE / SAFMIN
00429       CALL SLABAD( SAFMIN, SAFMAX )
00430       ULPINV = ONE / ULP
00431 *
00432 *     The values RMAGN(2:3) depend on N, see below.
00433 *
00434       RMAGN( 0 ) = ZERO
00435       RMAGN( 1 ) = ONE
00436 *
00437 *     Loop over matrix sizes
00438 *
00439       NTESTT = 0
00440       NERRS = 0
00441       NMATS = 0
00442 *
00443       DO 190 JSIZE = 1, NSIZES
00444          N = NN( JSIZE )
00445          N1 = MAX( 1, N )
00446          RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 )
00447          RMAGN( 3 ) = SAFMIN*ULPINV*REAL( N1 )
00448 *
00449          IF( NSIZES.NE.1 ) THEN
00450             MTYPES = MIN( MAXTYP, NTYPES )
00451          ELSE
00452             MTYPES = MIN( MAXTYP+1, NTYPES )
00453          END IF
00454 *
00455 *        Loop over matrix types
00456 *
00457          DO 180 JTYPE = 1, MTYPES
00458             IF( .NOT.DOTYPE( JTYPE ) )
00459      $         GO TO 180
00460             NMATS = NMATS + 1
00461             NTEST = 0
00462 *
00463 *           Save ISEED in case of an error.
00464 *
00465             DO 20 J = 1, 4
00466                IOLDSD( J ) = ISEED( J )
00467    20       CONTINUE
00468 *
00469 *           Initialize RESULT
00470 *
00471             DO 30 J = 1, 13
00472                RESULT( J ) = ZERO
00473    30       CONTINUE
00474 *
00475 *           Generate test matrices A and B
00476 *
00477 *           Description of control parameters:
00478 *
00479 *           KCLASS: =1 means w/o rotation, =2 means w/ rotation,
00480 *                   =3 means random.
00481 *           KATYPE: the "type" to be passed to SLATM4 for computing A.
00482 *           KAZERO: the pattern of zeros on the diagonal for A:
00483 *                   =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
00484 *                   =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
00485 *                   =6: ( 0, 1, 0, xxx, 0 ).  (xxx means a string of
00486 *                   non-zero entries.)
00487 *           KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
00488 *                   =2: large, =3: small.
00489 *           IASIGN: 1 if the diagonal elements of A are to be
00490 *                   multiplied by a random magnitude 1 number, =2 if
00491 *                   randomly chosen diagonal blocks are to be rotated
00492 *                   to form 2x2 blocks.
00493 *           KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
00494 *           KTRIAN: =0: don't fill in the upper triangle, =1: do.
00495 *           KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
00496 *           RMAGN: used to implement KAMAGN and KBMAGN.
00497 *
00498             IF( MTYPES.GT.MAXTYP )
00499      $         GO TO 110
00500             IINFO = 0
00501             IF( KCLASS( JTYPE ).LT.3 ) THEN
00502 *
00503 *              Generate A (w/o rotation)
00504 *
00505                IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
00506                   IN = 2*( ( N-1 ) / 2 ) + 1
00507                   IF( IN.NE.N )
00508      $               CALL SLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
00509                ELSE
00510                   IN = N
00511                END IF
00512                CALL SLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
00513      $                      KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ),
00514      $                      RMAGN( KAMAGN( JTYPE ) ), ULP,
00515      $                      RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
00516      $                      ISEED, A, LDA )
00517                IADD = KADD( KAZERO( JTYPE ) )
00518                IF( IADD.GT.0 .AND. IADD.LE.N )
00519      $            A( IADD, IADD ) = ONE
00520 *
00521 *              Generate B (w/o rotation)
00522 *
00523                IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
00524                   IN = 2*( ( N-1 ) / 2 ) + 1
00525                   IF( IN.NE.N )
00526      $               CALL SLASET( 'Full', N, N, ZERO, ZERO, B, LDA )
00527                ELSE
00528                   IN = N
00529                END IF
00530                CALL SLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
00531      $                      KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ),
00532      $                      RMAGN( KBMAGN( JTYPE ) ), ONE,
00533      $                      RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
00534      $                      ISEED, B, LDA )
00535                IADD = KADD( KBZERO( JTYPE ) )
00536                IF( IADD.NE.0 .AND. IADD.LE.N )
00537      $            B( IADD, IADD ) = ONE
00538 *
00539                IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
00540 *
00541 *                 Include rotations
00542 *
00543 *                 Generate Q, Z as Householder transformations times
00544 *                 a diagonal matrix.
00545 *
00546                   DO 50 JC = 1, N - 1
00547                      DO 40 JR = JC, N
00548                         Q( JR, JC ) = SLARND( 3, ISEED )
00549                         Z( JR, JC ) = SLARND( 3, ISEED )
00550    40                CONTINUE
00551                      CALL SLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
00552      $                            WORK( JC ) )
00553                      WORK( 2*N+JC ) = SIGN( ONE, Q( JC, JC ) )
00554                      Q( JC, JC ) = ONE
00555                      CALL SLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
00556      $                            WORK( N+JC ) )
00557                      WORK( 3*N+JC ) = SIGN( ONE, Z( JC, JC ) )
00558                      Z( JC, JC ) = ONE
00559    50             CONTINUE
00560                   Q( N, N ) = ONE
00561                   WORK( N ) = ZERO
00562                   WORK( 3*N ) = SIGN( ONE, SLARND( 2, ISEED ) )
00563                   Z( N, N ) = ONE
00564                   WORK( 2*N ) = ZERO
00565                   WORK( 4*N ) = SIGN( ONE, SLARND( 2, ISEED ) )
00566 *
00567 *                 Apply the diagonal matrices
00568 *
00569                   DO 70 JC = 1, N
00570                      DO 60 JR = 1, N
00571                         A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
00572      $                                A( JR, JC )
00573                         B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
00574      $                                B( JR, JC )
00575    60                CONTINUE
00576    70             CONTINUE
00577                   CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
00578      $                         LDA, WORK( 2*N+1 ), IINFO )
00579                   IF( IINFO.NE.0 )
00580      $               GO TO 100
00581                   CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
00582      $                         A, LDA, WORK( 2*N+1 ), IINFO )
00583                   IF( IINFO.NE.0 )
00584      $               GO TO 100
00585                   CALL SORM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
00586      $                         LDA, WORK( 2*N+1 ), IINFO )
00587                   IF( IINFO.NE.0 )
00588      $               GO TO 100
00589                   CALL SORM2R( 'R', 'T', N, N, N-1, Z, LDQ, WORK( N+1 ),
00590      $                         B, LDA, WORK( 2*N+1 ), IINFO )
00591                   IF( IINFO.NE.0 )
00592      $               GO TO 100
00593                END IF
00594             ELSE
00595 *
00596 *              Random matrices
00597 *
00598                DO 90 JC = 1, N
00599                   DO 80 JR = 1, N
00600                      A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
00601      $                             SLARND( 2, ISEED )
00602                      B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
00603      $                             SLARND( 2, ISEED )
00604    80             CONTINUE
00605    90          CONTINUE
00606             END IF
00607 *
00608   100       CONTINUE
00609 *
00610             IF( IINFO.NE.0 ) THEN
00611                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
00612      $            IOLDSD
00613                INFO = ABS( IINFO )
00614                RETURN
00615             END IF
00616 *
00617   110       CONTINUE
00618 *
00619             DO 120 I = 1, 13
00620                RESULT( I ) = -ONE
00621   120       CONTINUE
00622 *
00623 *           Test with and without sorting of eigenvalues
00624 *
00625             DO 150 ISORT = 0, 1
00626                IF( ISORT.EQ.0 ) THEN
00627                   SORT = 'N'
00628                   RSUB = 0
00629                ELSE
00630                   SORT = 'S'
00631                   RSUB = 5
00632                END IF
00633 *
00634 *              Call SGGES to compute H, T, Q, Z, alpha, and beta.
00635 *
00636                CALL SLACPY( 'Full', N, N, A, LDA, S, LDA )
00637                CALL SLACPY( 'Full', N, N, B, LDA, T, LDA )
00638                NTEST = 1 + RSUB + ISORT
00639                RESULT( 1+RSUB+ISORT ) = ULPINV
00640                CALL SGGES( 'V', 'V', SORT, SLCTES, N, S, LDA, T, LDA,
00641      $                     SDIM, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDQ,
00642      $                     WORK, LWORK, BWORK, IINFO )
00643                IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
00644                   RESULT( 1+RSUB+ISORT ) = ULPINV
00645                   WRITE( NOUNIT, FMT = 9999 )'SGGES', IINFO, N, JTYPE,
00646      $               IOLDSD
00647                   INFO = ABS( IINFO )
00648                   GO TO 160
00649                END IF
00650 *
00651                NTEST = 4 + RSUB
00652 *
00653 *              Do tests 1--4 (or tests 7--9 when reordering )
00654 *
00655                IF( ISORT.EQ.0 ) THEN
00656                   CALL SGET51( 1, N, A, LDA, S, LDA, Q, LDQ, Z, LDQ,
00657      $                         WORK, RESULT( 1 ) )
00658                   CALL SGET51( 1, N, B, LDA, T, LDA, Q, LDQ, Z, LDQ,
00659      $                         WORK, RESULT( 2 ) )
00660                ELSE
00661                   CALL SGET54( N, A, LDA, B, LDA, S, LDA, T, LDA, Q,
00662      $                         LDQ, Z, LDQ, WORK, RESULT( 7 ) )
00663                END IF
00664                CALL SGET51( 3, N, A, LDA, T, LDA, Q, LDQ, Q, LDQ, WORK,
00665      $                      RESULT( 3+RSUB ) )
00666                CALL SGET51( 3, N, B, LDA, T, LDA, Z, LDQ, Z, LDQ, WORK,
00667      $                      RESULT( 4+RSUB ) )
00668 *
00669 *              Do test 5 and 6 (or Tests 10 and 11 when reordering):
00670 *              check Schur form of A and compare eigenvalues with
00671 *              diagonals.
00672 *
00673                NTEST = 6 + RSUB
00674                TEMP1 = ZERO
00675 *
00676                DO 130 J = 1, N
00677                   ILABAD = .FALSE.
00678                   IF( ALPHAI( J ).EQ.ZERO ) THEN
00679                      TEMP2 = ( ABS( ALPHAR( J )-S( J, J ) ) /
00680      $                       MAX( SAFMIN, ABS( ALPHAR( J ) ), ABS( S( J,
00681      $                       J ) ) )+ABS( BETA( J )-T( J, J ) ) /
00682      $                       MAX( SAFMIN, ABS( BETA( J ) ), ABS( T( J,
00683      $                       J ) ) ) ) / ULP
00684 *
00685                      IF( J.LT.N ) THEN
00686                         IF( S( J+1, J ).NE.ZERO ) THEN
00687                            ILABAD = .TRUE.
00688                            RESULT( 5+RSUB ) = ULPINV
00689                         END IF
00690                      END IF
00691                      IF( J.GT.1 ) THEN
00692                         IF( S( J, J-1 ).NE.ZERO ) THEN
00693                            ILABAD = .TRUE.
00694                            RESULT( 5+RSUB ) = ULPINV
00695                         END IF
00696                      END IF
00697 *
00698                   ELSE
00699                      IF( ALPHAI( J ).GT.ZERO ) THEN
00700                         I1 = J
00701                      ELSE
00702                         I1 = J - 1
00703                      END IF
00704                      IF( I1.LE.0 .OR. I1.GE.N ) THEN
00705                         ILABAD = .TRUE.
00706                      ELSE IF( I1.LT.N-1 ) THEN
00707                         IF( S( I1+2, I1+1 ).NE.ZERO ) THEN
00708                            ILABAD = .TRUE.
00709                            RESULT( 5+RSUB ) = ULPINV
00710                         END IF
00711                      ELSE IF( I1.GT.1 ) THEN
00712                         IF( S( I1, I1-1 ).NE.ZERO ) THEN
00713                            ILABAD = .TRUE.
00714                            RESULT( 5+RSUB ) = ULPINV
00715                         END IF
00716                      END IF
00717                      IF( .NOT.ILABAD ) THEN
00718                         CALL SGET53( S( I1, I1 ), LDA, T( I1, I1 ), LDA,
00719      $                               BETA( J ), ALPHAR( J ),
00720      $                               ALPHAI( J ), TEMP2, IERR )
00721                         IF( IERR.GE.3 ) THEN
00722                            WRITE( NOUNIT, FMT = 9998 )IERR, J, N,
00723      $                        JTYPE, IOLDSD
00724                            INFO = ABS( IERR )
00725                         END IF
00726                      ELSE
00727                         TEMP2 = ULPINV
00728                      END IF
00729 *
00730                   END IF
00731                   TEMP1 = MAX( TEMP1, TEMP2 )
00732                   IF( ILABAD ) THEN
00733                      WRITE( NOUNIT, FMT = 9997 )J, N, JTYPE, IOLDSD
00734                   END IF
00735   130          CONTINUE
00736                RESULT( 6+RSUB ) = TEMP1
00737 *
00738                IF( ISORT.GE.1 ) THEN
00739 *
00740 *                 Do test 12
00741 *
00742                   NTEST = 12
00743                   RESULT( 12 ) = ZERO
00744                   KNTEIG = 0
00745                   DO 140 I = 1, N
00746                      IF( SLCTES( ALPHAR( I ), ALPHAI( I ),
00747      $                   BETA( I ) ) .OR. SLCTES( ALPHAR( I ),
00748      $                   -ALPHAI( I ), BETA( I ) ) ) THEN
00749                         KNTEIG = KNTEIG + 1
00750                      END IF
00751                      IF( I.LT.N ) THEN
00752                         IF( ( SLCTES( ALPHAR( I+1 ), ALPHAI( I+1 ),
00753      $                      BETA( I+1 ) ) .OR. SLCTES( ALPHAR( I+1 ),
00754      $                      -ALPHAI( I+1 ), BETA( I+1 ) ) ) .AND.
00755      $                      ( .NOT.( SLCTES( ALPHAR( I ), ALPHAI( I ),
00756      $                      BETA( I ) ) .OR. SLCTES( ALPHAR( I ),
00757      $                      -ALPHAI( I ), BETA( I ) ) ) ) .AND.
00758      $                      IINFO.NE.N+2 ) THEN
00759                            RESULT( 12 ) = ULPINV
00760                         END IF
00761                      END IF
00762   140             CONTINUE
00763                   IF( SDIM.NE.KNTEIG ) THEN
00764                      RESULT( 12 ) = ULPINV
00765                   END IF
00766                END IF
00767 *
00768   150       CONTINUE
00769 *
00770 *           End of Loop -- Check for RESULT(j) > THRESH
00771 *
00772   160       CONTINUE
00773 *
00774             NTESTT = NTESTT + NTEST
00775 *
00776 *           Print out tests which fail.
00777 *
00778             DO 170 JR = 1, NTEST
00779                IF( RESULT( JR ).GE.THRESH ) THEN
00780 *
00781 *                 If this is the first test to fail,
00782 *                 print a header to the data file.
00783 *
00784                   IF( NERRS.EQ.0 ) THEN
00785                      WRITE( NOUNIT, FMT = 9996 )'SGS'
00786 *
00787 *                    Matrix types
00788 *
00789                      WRITE( NOUNIT, FMT = 9995 )
00790                      WRITE( NOUNIT, FMT = 9994 )
00791                      WRITE( NOUNIT, FMT = 9993 )'Orthogonal'
00792 *
00793 *                    Tests performed
00794 *
00795                      WRITE( NOUNIT, FMT = 9992 )'orthogonal', '''',
00796      $                  'transpose', ( '''', J = 1, 8 )
00797 *
00798                   END IF
00799                   NERRS = NERRS + 1
00800                   IF( RESULT( JR ).LT.10000.0 ) THEN
00801                      WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
00802      $                  RESULT( JR )
00803                   ELSE
00804                      WRITE( NOUNIT, FMT = 9990 )N, JTYPE, IOLDSD, JR,
00805      $                  RESULT( JR )
00806                   END IF
00807                END IF
00808   170       CONTINUE
00809 *
00810   180    CONTINUE
00811   190 CONTINUE
00812 *
00813 *     Summary
00814 *
00815       CALL ALASVM( 'SGS', NOUNIT, NERRS, NTESTT, 0 )
00816 *
00817       WORK( 1 ) = MAXWRK
00818 *
00819       RETURN
00820 *
00821  9999 FORMAT( ' SDRGES: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00822      $      I6, ', JTYPE=', I6, ', ISEED=(', 4( I4, ',' ), I5, ')' )
00823 *
00824  9998 FORMAT( ' SDRGES: SGET53 returned INFO=', I1, ' for eigenvalue ',
00825      $      I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(',
00826      $      4( I4, ',' ), I5, ')' )
00827 *
00828  9997 FORMAT( ' SDRGES: S not in Schur form at eigenvalue ', I6, '.',
00829      $      / 9X, 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ),
00830      $      I5, ')' )
00831 *
00832  9996 FORMAT( / 1X, A3, ' -- Real Generalized Schur form driver' )
00833 *
00834  9995 FORMAT( ' Matrix types (see SDRGES for details): ' )
00835 *
00836  9994 FORMAT( ' Special Matrices:', 23X,
00837      $      '(J''=transposed Jordan block)',
00838      $      / '   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I)  5=(J'',J'')  ',
00839      $      '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices:  ( ',
00840      $      'D=diag(0,1,2,...) )', / '   7=(D,I)   9=(large*D, small*I',
00841      $      ')  11=(large*I, small*D)  13=(large*D, large*I)', /
00842      $      '   8=(I,D)  10=(small*D, large*I)  12=(small*I, large*D) ',
00843      $      ' 14=(small*D, small*I)', / '  15=(D, reversed D)' )
00844  9993 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
00845      $      / '  16=Transposed Jordan Blocks             19=geometric ',
00846      $      'alpha, beta=0,1', / '  17=arithm. alpha&beta             ',
00847      $      '      20=arithmetic alpha, beta=0,1', / '  18=clustered ',
00848      $      'alpha, beta=0,1            21=random alpha, beta=0,1',
00849      $      / ' Large & Small Matrices:', / '  22=(large, small)   ',
00850      $      '23=(small,large)    24=(small,small)    25=(large,large)',
00851      $      / '  26=random O(1) matrices.' )
00852 *
00853  9992 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00854      $      'Q and Z are ', A, ',', / 19X,
00855      $      'l and r are the appropriate left and right', / 19X,
00856      $      'eigenvectors, resp., a is alpha, b is beta, and', / 19X, A,
00857      $      ' means ', A, '.)', / ' Without ordering: ',
00858      $      / '  1 = | A - Q S Z', A,
00859      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
00860      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
00861      $      ' | / ( n ulp )             4 = | I - ZZ', A,
00862      $      ' | / ( n ulp )', / '  5 = A is in Schur form S',
00863      $      / '  6 = difference between (alpha,beta)',
00864      $      ' and diagonals of (S,T)', / ' With ordering: ',
00865      $      / '  7 = | (A,B) - Q (S,T) Z', A,
00866      $      ' | / ( |(A,B)| n ulp )  ', / '  8 = | I - QQ', A,
00867      $      ' | / ( n ulp )            9 = | I - ZZ', A,
00868      $      ' | / ( n ulp )', / ' 10 = A is in Schur form S',
00869      $      / ' 11 = difference between (alpha,beta) and diagonals',
00870      $      ' of (S,T)', / ' 12 = SDIM is the correct number of ',
00871      $      'selected eigenvalues', / )
00872  9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
00873      $      4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
00874  9990 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
00875      $      4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
00876 *
00877 *     End of SDRGES
00878 *
00879       END