LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
schkgg.f
Go to the documentation of this file.
00001 *> \brief \b SCHKGG
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE SCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00012 *                          TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
00013 *                          S2, P1, P2, U, LDU, V, Q, Z, ALPHR1, ALPHI1,
00014 *                          BETA1, ALPHR3, ALPHI3, BETA3, EVECTL, EVECTR,
00015 *                          WORK, LWORK, LLWORK, RESULT, INFO )
00016 * 
00017 *       .. Scalar Arguments ..
00018 *       LOGICAL            TSTDIF
00019 *       INTEGER            INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
00020 *       REAL               THRESH, THRSHN
00021 *       ..
00022 *       .. Array Arguments ..
00023 *       LOGICAL            DOTYPE( * ), LLWORK( * )
00024 *       INTEGER            ISEED( 4 ), NN( * )
00025 *       REAL               A( LDA, * ), ALPHI1( * ), ALPHI3( * ),
00026 *      $                   ALPHR1( * ), ALPHR3( * ), B( LDA, * ),
00027 *      $                   BETA1( * ), BETA3( * ), EVECTL( LDU, * ),
00028 *      $                   EVECTR( LDU, * ), H( LDA, * ), P1( LDA, * ),
00029 *      $                   P2( LDA, * ), Q( LDU, * ), RESULT( 15 ),
00030 *      $                   S1( LDA, * ), S2( LDA, * ), T( LDA, * ),
00031 *      $                   U( LDU, * ), V( LDU, * ), WORK( * ),
00032 *      $                   Z( LDU, * )
00033 *       ..
00034 *  
00035 *
00036 *> \par Purpose:
00037 *  =============
00038 *>
00039 *> \verbatim
00040 *>
00041 *> SCHKGG  checks the nonsymmetric generalized eigenvalue problem
00042 *> routines.
00043 *>                                T          T        T
00044 *> SGGHRD factors A and B as U H V  and U T V , where   means
00045 *> transpose, H is hessenberg, T is triangular and U and V are
00046 *> orthogonal.
00047 *>                                 T          T
00048 *> SHGEQZ factors H and T as  Q S Z  and Q P Z , where P is upper
00049 *> triangular, S is in generalized Schur form (block upper triangular,
00050 *> with 1x1 and 2x2 blocks on the diagonal, the 2x2 blocks
00051 *> corresponding to complex conjugate pairs of generalized
00052 *> eigenvalues), and Q and Z are orthogonal.  It also computes the
00053 *> generalized eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)),
00054 *> where alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus,
00055 *> w(j) = alpha(j)/beta(j) is a root of the generalized eigenvalue
00056 *> problem
00057 *>
00058 *>     det( A - w(j) B ) = 0
00059 *>
00060 *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
00061 *> problem
00062 *>
00063 *>     det( m(j) A - B ) = 0
00064 *>
00065 *> STGEVC computes the matrix L of left eigenvectors and the matrix R
00066 *> of right eigenvectors for the matrix pair ( S, P ).  In the
00067 *> description below,  l and r are left and right eigenvectors
00068 *> corresponding to the generalized eigenvalues (alpha,beta).
00069 *>
00070 *> When SCHKGG is called, a number of matrix "sizes" ("n's") and a
00071 *> number of matrix "types" are specified.  For each size ("n")
00072 *> and each type of matrix, one matrix will be generated and used
00073 *> to test the nonsymmetric eigenroutines.  For each matrix, 15
00074 *> tests will be performed.  The first twelve "test ratios" should be
00075 *> small -- O(1).  They will be compared with the threshhold THRESH:
00076 *>
00077 *>                  T
00078 *> (1)   | A - U H V  | / ( |A| n ulp )
00079 *>
00080 *>                  T
00081 *> (2)   | B - U T V  | / ( |B| n ulp )
00082 *>
00083 *>               T
00084 *> (3)   | I - UU  | / ( n ulp )
00085 *>
00086 *>               T
00087 *> (4)   | I - VV  | / ( n ulp )
00088 *>
00089 *>                  T
00090 *> (5)   | H - Q S Z  | / ( |H| n ulp )
00091 *>
00092 *>                  T
00093 *> (6)   | T - Q P Z  | / ( |T| n ulp )
00094 *>
00095 *>               T
00096 *> (7)   | I - QQ  | / ( n ulp )
00097 *>
00098 *>               T
00099 *> (8)   | I - ZZ  | / ( n ulp )
00100 *>
00101 *> (9)   max over all left eigenvalue/-vector pairs (beta/alpha,l) of
00102 *>
00103 *>    | l**H * (beta S - alpha P) | / ( ulp max( |beta S|, |alpha P| ) )
00104 *>
00105 *> (10)  max over all left eigenvalue/-vector pairs (beta/alpha,l') of
00106 *>                           T
00107 *>   | l'**H * (beta H - alpha T) | / ( ulp max( |beta H|, |alpha T| ) )
00108 *>
00109 *>       where the eigenvectors l' are the result of passing Q to
00110 *>       STGEVC and back transforming (HOWMNY='B').
00111 *>
00112 *> (11)  max over all right eigenvalue/-vector pairs (beta/alpha,r) of
00113 *>
00114 *>       | (beta S - alpha T) r | / ( ulp max( |beta S|, |alpha T| ) )
00115 *>
00116 *> (12)  max over all right eigenvalue/-vector pairs (beta/alpha,r') of
00117 *>
00118 *>       | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
00119 *>
00120 *>       where the eigenvectors r' are the result of passing Z to
00121 *>       STGEVC and back transforming (HOWMNY='B').
00122 *>
00123 *> The last three test ratios will usually be small, but there is no
00124 *> mathematical requirement that they be so.  They are therefore
00125 *> compared with THRESH only if TSTDIF is .TRUE.
00126 *>
00127 *> (13)  | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
00128 *>
00129 *> (14)  | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
00130 *>
00131 *> (15)  max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
00132 *>            |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
00133 *>
00134 *> In addition, the normalization of L and R are checked, and compared
00135 *> with the threshhold THRSHN.
00136 *>
00137 *> Test Matrices
00138 *> ---- --------
00139 *>
00140 *> The sizes of the test matrices are specified by an array
00141 *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
00142 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
00143 *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
00144 *> Currently, the list of possible types is:
00145 *>
00146 *> (1)  ( 0, 0 )         (a pair of zero matrices)
00147 *>
00148 *> (2)  ( I, 0 )         (an identity and a zero matrix)
00149 *>
00150 *> (3)  ( 0, I )         (an identity and a zero matrix)
00151 *>
00152 *> (4)  ( I, I )         (a pair of identity matrices)
00153 *>
00154 *>         t   t
00155 *> (5)  ( J , J  )       (a pair of transposed Jordan blocks)
00156 *>
00157 *>                                     t                ( I   0  )
00158 *> (6)  ( X, Y )         where  X = ( J   0  )  and Y = (      t )
00159 *>                                  ( 0   I  )          ( 0   J  )
00160 *>                       and I is a k x k identity and J a (k+1)x(k+1)
00161 *>                       Jordan block; k=(N-1)/2
00162 *>
00163 *> (7)  ( D, I )         where D is diag( 0, 1,..., N-1 ) (a diagonal
00164 *>                       matrix with those diagonal entries.)
00165 *> (8)  ( I, D )
00166 *>
00167 *> (9)  ( big*D, small*I ) where "big" is near overflow and small=1/big
00168 *>
00169 *> (10) ( small*D, big*I )
00170 *>
00171 *> (11) ( big*I, small*D )
00172 *>
00173 *> (12) ( small*I, big*D )
00174 *>
00175 *> (13) ( big*D, big*I )
00176 *>
00177 *> (14) ( small*D, small*I )
00178 *>
00179 *> (15) ( D1, D2 )        where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
00180 *>                        D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
00181 *>           t   t
00182 *> (16) U ( J , J ) V     where U and V are random orthogonal matrices.
00183 *>
00184 *> (17) U ( T1, T2 ) V    where T1 and T2 are upper triangular matrices
00185 *>                        with random O(1) entries above the diagonal
00186 *>                        and diagonal entries diag(T1) =
00187 *>                        ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
00188 *>                        ( 0, N-3, N-4,..., 1, 0, 0 )
00189 *>
00190 *> (18) U ( T1, T2 ) V    diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
00191 *>                        diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
00192 *>                        s = machine precision.
00193 *>
00194 *> (19) U ( T1, T2 ) V    diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
00195 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
00196 *>
00197 *>                                                        N-5
00198 *> (20) U ( T1, T2 ) V    diag(T1)=( 0, 0, 1, 1, a, ..., a   =s, 0 )
00199 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
00200 *>
00201 *> (21) U ( T1, T2 ) V    diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
00202 *>                        diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
00203 *>                        where r1,..., r(N-4) are random.
00204 *>
00205 *> (22) U ( big*T1, small*T2 ) V    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00206 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00207 *>
00208 *> (23) U ( small*T1, big*T2 ) V    diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00209 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00210 *>
00211 *> (24) U ( small*T1, small*T2 ) V  diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00212 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00213 *>
00214 *> (25) U ( big*T1, big*T2 ) V      diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
00215 *>                                  diag(T2) = ( 0, 1, ..., 1, 0, 0 )
00216 *>
00217 *> (26) U ( T1, T2 ) V     where T1 and T2 are random upper-triangular
00218 *>                         matrices.
00219 *> \endverbatim
00220 *
00221 *  Arguments:
00222 *  ==========
00223 *
00224 *> \param[in] NSIZES
00225 *> \verbatim
00226 *>          NSIZES is INTEGER
00227 *>          The number of sizes of matrices to use.  If it is zero,
00228 *>          SCHKGG does nothing.  It must be at least zero.
00229 *> \endverbatim
00230 *>
00231 *> \param[in] NN
00232 *> \verbatim
00233 *>          NN is INTEGER array, dimension (NSIZES)
00234 *>          An array containing the sizes to be used for the matrices.
00235 *>          Zero values will be skipped.  The values must be at least
00236 *>          zero.
00237 *> \endverbatim
00238 *>
00239 *> \param[in] NTYPES
00240 *> \verbatim
00241 *>          NTYPES is INTEGER
00242 *>          The number of elements in DOTYPE.   If it is zero, SCHKGG
00243 *>          does nothing.  It must be at least zero.  If it is MAXTYP+1
00244 *>          and NSIZES is 1, then an additional type, MAXTYP+1 is
00245 *>          defined, which is to use whatever matrix is in A.  This
00246 *>          is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
00247 *>          DOTYPE(MAXTYP+1) is .TRUE. .
00248 *> \endverbatim
00249 *>
00250 *> \param[in] DOTYPE
00251 *> \verbatim
00252 *>          DOTYPE is LOGICAL array, dimension (NTYPES)
00253 *>          If DOTYPE(j) is .TRUE., then for each size in NN a
00254 *>          matrix of that size and of type j will be generated.
00255 *>          If NTYPES is smaller than the maximum number of types
00256 *>          defined (PARAMETER MAXTYP), then types NTYPES+1 through
00257 *>          MAXTYP will not be generated.  If NTYPES is larger
00258 *>          than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
00259 *>          will be ignored.
00260 *> \endverbatim
00261 *>
00262 *> \param[in,out] ISEED
00263 *> \verbatim
00264 *>          ISEED is INTEGER array, dimension (4)
00265 *>          On entry ISEED specifies the seed of the random number
00266 *>          generator. The array elements should be between 0 and 4095;
00267 *>          if not they will be reduced mod 4096.  Also, ISEED(4) must
00268 *>          be odd.  The random number generator uses a linear
00269 *>          congruential sequence limited to small integers, and so
00270 *>          should produce machine independent random numbers. The
00271 *>          values of ISEED are changed on exit, and can be used in the
00272 *>          next call to SCHKGG to continue the same random number
00273 *>          sequence.
00274 *> \endverbatim
00275 *>
00276 *> \param[in] THRESH
00277 *> \verbatim
00278 *>          THRESH is REAL
00279 *>          A test will count as "failed" if the "error", computed as
00280 *>          described above, exceeds THRESH.  Note that the error is
00281 *>          scaled to be O(1), so THRESH should be a reasonably small
00282 *>          multiple of 1, e.g., 10 or 100.  In particular, it should
00283 *>          not depend on the precision (single vs. double) or the size
00284 *>          of the matrix.  It must be at least zero.
00285 *> \endverbatim
00286 *>
00287 *> \param[in] TSTDIF
00288 *> \verbatim
00289 *>          TSTDIF is LOGICAL
00290 *>          Specifies whether test ratios 13-15 will be computed and
00291 *>          compared with THRESH.
00292 *>          = .FALSE.: Only test ratios 1-12 will be computed and tested.
00293 *>                     Ratios 13-15 will be set to zero.
00294 *>          = .TRUE.:  All the test ratios 1-15 will be computed and
00295 *>                     tested.
00296 *> \endverbatim
00297 *>
00298 *> \param[in] THRSHN
00299 *> \verbatim
00300 *>          THRSHN is REAL
00301 *>          Threshhold for reporting eigenvector normalization error.
00302 *>          If the normalization of any eigenvector differs from 1 by
00303 *>          more than THRSHN*ulp, then a special error message will be
00304 *>          printed.  (This is handled separately from the other tests,
00305 *>          since only a compiler or programming error should cause an
00306 *>          error message, at least if THRSHN is at least 5--10.)
00307 *> \endverbatim
00308 *>
00309 *> \param[in] NOUNIT
00310 *> \verbatim
00311 *>          NOUNIT is INTEGER
00312 *>          The FORTRAN unit number for printing out error messages
00313 *>          (e.g., if a routine returns IINFO not equal to 0.)
00314 *> \endverbatim
00315 *>
00316 *> \param[in,out] A
00317 *> \verbatim
00318 *>          A is REAL array, dimension
00319 *>                            (LDA, max(NN))
00320 *>          Used to hold the original A matrix.  Used as input only
00321 *>          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
00322 *>          DOTYPE(MAXTYP+1)=.TRUE.
00323 *> \endverbatim
00324 *>
00325 *> \param[in] LDA
00326 *> \verbatim
00327 *>          LDA is INTEGER
00328 *>          The leading dimension of A, B, H, T, S1, P1, S2, and P2.
00329 *>          It must be at least 1 and at least max( NN ).
00330 *> \endverbatim
00331 *>
00332 *> \param[in,out] B
00333 *> \verbatim
00334 *>          B is REAL array, dimension
00335 *>                            (LDA, max(NN))
00336 *>          Used to hold the original B matrix.  Used as input only
00337 *>          if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
00338 *>          DOTYPE(MAXTYP+1)=.TRUE.
00339 *> \endverbatim
00340 *>
00341 *> \param[out] H
00342 *> \verbatim
00343 *>          H is REAL array, dimension (LDA, max(NN))
00344 *>          The upper Hessenberg matrix computed from A by SGGHRD.
00345 *> \endverbatim
00346 *>
00347 *> \param[out] T
00348 *> \verbatim
00349 *>          T is REAL array, dimension (LDA, max(NN))
00350 *>          The upper triangular matrix computed from B by SGGHRD.
00351 *> \endverbatim
00352 *>
00353 *> \param[out] S1
00354 *> \verbatim
00355 *>          S1 is REAL array, dimension (LDA, max(NN))
00356 *>          The Schur (block upper triangular) matrix computed from H by
00357 *>          SHGEQZ when Q and Z are also computed.
00358 *> \endverbatim
00359 *>
00360 *> \param[out] S2
00361 *> \verbatim
00362 *>          S2 is REAL array, dimension (LDA, max(NN))
00363 *>          The Schur (block upper triangular) matrix computed from H by
00364 *>          SHGEQZ when Q and Z are not computed.
00365 *> \endverbatim
00366 *>
00367 *> \param[out] P1
00368 *> \verbatim
00369 *>          P1 is REAL array, dimension (LDA, max(NN))
00370 *>          The upper triangular matrix computed from T by SHGEQZ
00371 *>          when Q and Z are also computed.
00372 *> \endverbatim
00373 *>
00374 *> \param[out] P2
00375 *> \verbatim
00376 *>          P2 is REAL array, dimension (LDA, max(NN))
00377 *>          The upper triangular matrix computed from T by SHGEQZ
00378 *>          when Q and Z are not computed.
00379 *> \endverbatim
00380 *>
00381 *> \param[out] U
00382 *> \verbatim
00383 *>          U is REAL array, dimension (LDU, max(NN))
00384 *>          The (left) orthogonal matrix computed by SGGHRD.
00385 *> \endverbatim
00386 *>
00387 *> \param[in] LDU
00388 *> \verbatim
00389 *>          LDU is INTEGER
00390 *>          The leading dimension of U, V, Q, Z, EVECTL, and EVECTR.  It
00391 *>          must be at least 1 and at least max( NN ).
00392 *> \endverbatim
00393 *>
00394 *> \param[out] V
00395 *> \verbatim
00396 *>          V is REAL array, dimension (LDU, max(NN))
00397 *>          The (right) orthogonal matrix computed by SGGHRD.
00398 *> \endverbatim
00399 *>
00400 *> \param[out] Q
00401 *> \verbatim
00402 *>          Q is REAL array, dimension (LDU, max(NN))
00403 *>          The (left) orthogonal matrix computed by SHGEQZ.
00404 *> \endverbatim
00405 *>
00406 *> \param[out] Z
00407 *> \verbatim
00408 *>          Z is REAL array, dimension (LDU, max(NN))
00409 *>          The (left) orthogonal matrix computed by SHGEQZ.
00410 *> \endverbatim
00411 *>
00412 *> \param[out] ALPHR1
00413 *> \verbatim
00414 *>          ALPHR1 is REAL array, dimension (max(NN))
00415 *> \endverbatim
00416 *>
00417 *> \param[out] ALPHI1
00418 *> \verbatim
00419 *>          ALPHI1 is REAL array, dimension (max(NN))
00420 *> \endverbatim
00421 *>
00422 *> \param[out] BETA1
00423 *> \verbatim
00424 *>          BETA1 is REAL array, dimension (max(NN))
00425 *>
00426 *>          The generalized eigenvalues of (A,B) computed by SHGEQZ
00427 *>          when Q, Z, and the full Schur matrices are computed.
00428 *>          On exit, ( ALPHR1(k)+ALPHI1(k)*i ) / BETA1(k) is the k-th
00429 *>          generalized eigenvalue of the matrices in A and B.
00430 *> \endverbatim
00431 *>
00432 *> \param[out] ALPHR3
00433 *> \verbatim
00434 *>          ALPHR3 is REAL array, dimension (max(NN))
00435 *> \endverbatim
00436 *>
00437 *> \param[out] ALPHI3
00438 *> \verbatim
00439 *>          ALPHI3 is REAL array, dimension (max(NN))
00440 *> \endverbatim
00441 *>
00442 *> \param[out] BETA3
00443 *> \verbatim
00444 *>          BETA3 is REAL array, dimension (max(NN))
00445 *> \endverbatim
00446 *>
00447 *> \param[out] EVECTL
00448 *> \verbatim
00449 *>          EVECTL is REAL array, dimension (LDU, max(NN))
00450 *>          The (block lower triangular) left eigenvector matrix for
00451 *>          the matrices in S1 and P1.  (See STGEVC for the format.)
00452 *> \endverbatim
00453 *>
00454 *> \param[out] EVECTR
00455 *> \verbatim
00456 *>          EVECTR is REAL array, dimension (LDU, max(NN))
00457 *>          The (block upper triangular) right eigenvector matrix for
00458 *>          the matrices in S1 and P1.  (See STGEVC for the format.)
00459 *> \endverbatim
00460 *>
00461 *> \param[out] WORK
00462 *> \verbatim
00463 *>          WORK is REAL array, dimension (LWORK)
00464 *> \endverbatim
00465 *>
00466 *> \param[in] LWORK
00467 *> \verbatim
00468 *>          LWORK is INTEGER
00469 *>          The number of entries in WORK.  This must be at least
00470 *>          max( 2 * N**2, 6*N, 1 ), for all N=NN(j).
00471 *> \endverbatim
00472 *>
00473 *> \param[out] LLWORK
00474 *> \verbatim
00475 *>          LLWORK is LOGICAL array, dimension (max(NN))
00476 *> \endverbatim
00477 *>
00478 *> \param[out] RESULT
00479 *> \verbatim
00480 *>          RESULT is REAL array, dimension (15)
00481 *>          The values computed by the tests described above.
00482 *>          The values are currently limited to 1/ulp, to avoid
00483 *>          overflow.
00484 *> \endverbatim
00485 *>
00486 *> \param[out] INFO
00487 *> \verbatim
00488 *>          INFO is INTEGER
00489 *>          = 0:  successful exit
00490 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00491 *>          > 0:  A routine returned an error code.  INFO is the
00492 *>                absolute value of the INFO value returned.
00493 *> \endverbatim
00494 *
00495 *  Authors:
00496 *  ========
00497 *
00498 *> \author Univ. of Tennessee 
00499 *> \author Univ. of California Berkeley 
00500 *> \author Univ. of Colorado Denver 
00501 *> \author NAG Ltd. 
00502 *
00503 *> \date November 2011
00504 *
00505 *> \ingroup single_eig
00506 *
00507 *  =====================================================================
00508       SUBROUTINE SCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
00509      $                   TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
00510      $                   S2, P1, P2, U, LDU, V, Q, Z, ALPHR1, ALPHI1,
00511      $                   BETA1, ALPHR3, ALPHI3, BETA3, EVECTL, EVECTR,
00512      $                   WORK, LWORK, LLWORK, RESULT, INFO )
00513 *
00514 *  -- LAPACK test routine (version 3.4.0) --
00515 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00516 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00517 *     November 2011
00518 *
00519 *     .. Scalar Arguments ..
00520       LOGICAL            TSTDIF
00521       INTEGER            INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
00522       REAL               THRESH, THRSHN
00523 *     ..
00524 *     .. Array Arguments ..
00525       LOGICAL            DOTYPE( * ), LLWORK( * )
00526       INTEGER            ISEED( 4 ), NN( * )
00527       REAL               A( LDA, * ), ALPHI1( * ), ALPHI3( * ),
00528      $                   ALPHR1( * ), ALPHR3( * ), B( LDA, * ),
00529      $                   BETA1( * ), BETA3( * ), EVECTL( LDU, * ),
00530      $                   EVECTR( LDU, * ), H( LDA, * ), P1( LDA, * ),
00531      $                   P2( LDA, * ), Q( LDU, * ), RESULT( 15 ),
00532      $                   S1( LDA, * ), S2( LDA, * ), T( LDA, * ),
00533      $                   U( LDU, * ), V( LDU, * ), WORK( * ),
00534      $                   Z( LDU, * )
00535 *     ..
00536 *
00537 *  =====================================================================
00538 *
00539 *     .. Parameters ..
00540       REAL               ZERO, ONE
00541       PARAMETER          ( ZERO = 0.0, ONE = 1.0 )
00542       INTEGER            MAXTYP
00543       PARAMETER          ( MAXTYP = 26 )
00544 *     ..
00545 *     .. Local Scalars ..
00546       LOGICAL            BADNN
00547       INTEGER            I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
00548      $                   LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
00549      $                   NTEST, NTESTT
00550       REAL               ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
00551      $                   ULP, ULPINV
00552 *     ..
00553 *     .. Local Arrays ..
00554       INTEGER            IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
00555      $                   IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
00556      $                   KATYPE( MAXTYP ), KAZERO( MAXTYP ),
00557      $                   KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
00558      $                   KBZERO( MAXTYP ), KCLASS( MAXTYP ),
00559      $                   KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
00560       REAL               DUMMA( 4 ), RMAGN( 0: 3 )
00561 *     ..
00562 *     .. External Functions ..
00563       REAL               SLAMCH, SLANGE, SLARND
00564       EXTERNAL           SLAMCH, SLANGE, SLARND
00565 *     ..
00566 *     .. External Subroutines ..
00567       EXTERNAL           SGEQR2, SGET51, SGET52, SGGHRD, SHGEQZ, SLABAD,
00568      $                   SLACPY, SLARFG, SLASET, SLASUM, SLATM4, SORM2R,
00569      $                   STGEVC, XERBLA
00570 *     ..
00571 *     .. Intrinsic Functions ..
00572       INTRINSIC          ABS, MAX, MIN, REAL, SIGN
00573 *     ..
00574 *     .. Data statements ..
00575       DATA               KCLASS / 15*1, 10*2, 1*3 /
00576       DATA               KZ1 / 0, 1, 2, 1, 3, 3 /
00577       DATA               KZ2 / 0, 0, 1, 2, 1, 1 /
00578       DATA               KADD / 0, 0, 0, 0, 3, 2 /
00579       DATA               KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
00580      $                   4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
00581       DATA               KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
00582      $                   1, 1, -4, 2, -4, 8*8, 0 /
00583       DATA               KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
00584      $                   4*5, 4*3, 1 /
00585       DATA               KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
00586      $                   4*6, 4*4, 1 /
00587       DATA               KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
00588      $                   2, 1 /
00589       DATA               KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
00590      $                   2, 1 /
00591       DATA               KTRIAN / 16*0, 10*1 /
00592       DATA               IASIGN / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
00593      $                   5*2, 0 /
00594       DATA               IBSIGN / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
00595 *     ..
00596 *     .. Executable Statements ..
00597 *
00598 *     Check for errors
00599 *
00600       INFO = 0
00601 *
00602       BADNN = .FALSE.
00603       NMAX = 1
00604       DO 10 J = 1, NSIZES
00605          NMAX = MAX( NMAX, NN( J ) )
00606          IF( NN( J ).LT.0 )
00607      $      BADNN = .TRUE.
00608    10 CONTINUE
00609 *
00610 *     Maximum blocksize and shift -- we assume that blocksize and number
00611 *     of shifts are monotone increasing functions of N.
00612 *
00613       LWKOPT = MAX( 6*NMAX, 2*NMAX*NMAX, 1 )
00614 *
00615 *     Check for errors
00616 *
00617       IF( NSIZES.LT.0 ) THEN
00618          INFO = -1
00619       ELSE IF( BADNN ) THEN
00620          INFO = -2
00621       ELSE IF( NTYPES.LT.0 ) THEN
00622          INFO = -3
00623       ELSE IF( THRESH.LT.ZERO ) THEN
00624          INFO = -6
00625       ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
00626          INFO = -10
00627       ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
00628          INFO = -19
00629       ELSE IF( LWKOPT.GT.LWORK ) THEN
00630          INFO = -30
00631       END IF
00632 *
00633       IF( INFO.NE.0 ) THEN
00634          CALL XERBLA( 'SCHKGG', -INFO )
00635          RETURN
00636       END IF
00637 *
00638 *     Quick return if possible
00639 *
00640       IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
00641      $   RETURN
00642 *
00643       SAFMIN = SLAMCH( 'Safe minimum' )
00644       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00645       SAFMIN = SAFMIN / ULP
00646       SAFMAX = ONE / SAFMIN
00647       CALL SLABAD( SAFMIN, SAFMAX )
00648       ULPINV = ONE / ULP
00649 *
00650 *     The values RMAGN(2:3) depend on N, see below.
00651 *
00652       RMAGN( 0 ) = ZERO
00653       RMAGN( 1 ) = ONE
00654 *
00655 *     Loop over sizes, types
00656 *
00657       NTESTT = 0
00658       NERRS = 0
00659       NMATS = 0
00660 *
00661       DO 240 JSIZE = 1, NSIZES
00662          N = NN( JSIZE )
00663          N1 = MAX( 1, N )
00664          RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 )
00665          RMAGN( 3 ) = SAFMIN*ULPINV*N1
00666 *
00667          IF( NSIZES.NE.1 ) THEN
00668             MTYPES = MIN( MAXTYP, NTYPES )
00669          ELSE
00670             MTYPES = MIN( MAXTYP+1, NTYPES )
00671          END IF
00672 *
00673          DO 230 JTYPE = 1, MTYPES
00674             IF( .NOT.DOTYPE( JTYPE ) )
00675      $         GO TO 230
00676             NMATS = NMATS + 1
00677             NTEST = 0
00678 *
00679 *           Save ISEED in case of an error.
00680 *
00681             DO 20 J = 1, 4
00682                IOLDSD( J ) = ISEED( J )
00683    20       CONTINUE
00684 *
00685 *           Initialize RESULT
00686 *
00687             DO 30 J = 1, 15
00688                RESULT( J ) = ZERO
00689    30       CONTINUE
00690 *
00691 *           Compute A and B
00692 *
00693 *           Description of control parameters:
00694 *
00695 *           KCLASS: =1 means w/o rotation, =2 means w/ rotation,
00696 *                   =3 means random.
00697 *           KATYPE: the "type" to be passed to SLATM4 for computing A.
00698 *           KAZERO: the pattern of zeros on the diagonal for A:
00699 *                   =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
00700 *                   =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
00701 *                   =6: ( 0, 1, 0, xxx, 0 ).  (xxx means a string of
00702 *                   non-zero entries.)
00703 *           KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
00704 *                   =2: large, =3: small.
00705 *           IASIGN: 1 if the diagonal elements of A are to be
00706 *                   multiplied by a random magnitude 1 number, =2 if
00707 *                   randomly chosen diagonal blocks are to be rotated
00708 *                   to form 2x2 blocks.
00709 *           KBTYPE, KBZERO, KBMAGN, IBSIGN: the same, but for B.
00710 *           KTRIAN: =0: don't fill in the upper triangle, =1: do.
00711 *           KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
00712 *           RMAGN: used to implement KAMAGN and KBMAGN.
00713 *
00714             IF( MTYPES.GT.MAXTYP )
00715      $         GO TO 110
00716             IINFO = 0
00717             IF( KCLASS( JTYPE ).LT.3 ) THEN
00718 *
00719 *              Generate A (w/o rotation)
00720 *
00721                IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
00722                   IN = 2*( ( N-1 ) / 2 ) + 1
00723                   IF( IN.NE.N )
00724      $               CALL SLASET( 'Full', N, N, ZERO, ZERO, A, LDA )
00725                ELSE
00726                   IN = N
00727                END IF
00728                CALL SLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
00729      $                      KZ2( KAZERO( JTYPE ) ), IASIGN( JTYPE ),
00730      $                      RMAGN( KAMAGN( JTYPE ) ), ULP,
00731      $                      RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
00732      $                      ISEED, A, LDA )
00733                IADD = KADD( KAZERO( JTYPE ) )
00734                IF( IADD.GT.0 .AND. IADD.LE.N )
00735      $            A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
00736 *
00737 *              Generate B (w/o rotation)
00738 *
00739                IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
00740                   IN = 2*( ( N-1 ) / 2 ) + 1
00741                   IF( IN.NE.N )
00742      $               CALL SLASET( 'Full', N, N, ZERO, ZERO, B, LDA )
00743                ELSE
00744                   IN = N
00745                END IF
00746                CALL SLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
00747      $                      KZ2( KBZERO( JTYPE ) ), IBSIGN( JTYPE ),
00748      $                      RMAGN( KBMAGN( JTYPE ) ), ONE,
00749      $                      RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
00750      $                      ISEED, B, LDA )
00751                IADD = KADD( KBZERO( JTYPE ) )
00752                IF( IADD.NE.0 .AND. IADD.LE.N )
00753      $            B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
00754 *
00755                IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
00756 *
00757 *                 Include rotations
00758 *
00759 *                 Generate U, V as Householder transformations times
00760 *                 a diagonal matrix.
00761 *
00762                   DO 50 JC = 1, N - 1
00763                      DO 40 JR = JC, N
00764                         U( JR, JC ) = SLARND( 3, ISEED )
00765                         V( JR, JC ) = SLARND( 3, ISEED )
00766    40                CONTINUE
00767                      CALL SLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
00768      $                            WORK( JC ) )
00769                      WORK( 2*N+JC ) = SIGN( ONE, U( JC, JC ) )
00770                      U( JC, JC ) = ONE
00771                      CALL SLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
00772      $                            WORK( N+JC ) )
00773                      WORK( 3*N+JC ) = SIGN( ONE, V( JC, JC ) )
00774                      V( JC, JC ) = ONE
00775    50             CONTINUE
00776                   U( N, N ) = ONE
00777                   WORK( N ) = ZERO
00778                   WORK( 3*N ) = SIGN( ONE, SLARND( 2, ISEED ) )
00779                   V( N, N ) = ONE
00780                   WORK( 2*N ) = ZERO
00781                   WORK( 4*N ) = SIGN( ONE, SLARND( 2, ISEED ) )
00782 *
00783 *                 Apply the diagonal matrices
00784 *
00785                   DO 70 JC = 1, N
00786                      DO 60 JR = 1, N
00787                         A( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
00788      $                                A( JR, JC )
00789                         B( JR, JC ) = WORK( 2*N+JR )*WORK( 3*N+JC )*
00790      $                                B( JR, JC )
00791    60                CONTINUE
00792    70             CONTINUE
00793                   CALL SORM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A,
00794      $                         LDA, WORK( 2*N+1 ), IINFO )
00795                   IF( IINFO.NE.0 )
00796      $               GO TO 100
00797                   CALL SORM2R( 'R', 'T', N, N, N-1, V, LDU, WORK( N+1 ),
00798      $                         A, LDA, WORK( 2*N+1 ), IINFO )
00799                   IF( IINFO.NE.0 )
00800      $               GO TO 100
00801                   CALL SORM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B,
00802      $                         LDA, WORK( 2*N+1 ), IINFO )
00803                   IF( IINFO.NE.0 )
00804      $               GO TO 100
00805                   CALL SORM2R( 'R', 'T', N, N, N-1, V, LDU, WORK( N+1 ),
00806      $                         B, LDA, WORK( 2*N+1 ), IINFO )
00807                   IF( IINFO.NE.0 )
00808      $               GO TO 100
00809                END IF
00810             ELSE
00811 *
00812 *              Random matrices
00813 *
00814                DO 90 JC = 1, N
00815                   DO 80 JR = 1, N
00816                      A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
00817      $                             SLARND( 2, ISEED )
00818                      B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
00819      $                             SLARND( 2, ISEED )
00820    80             CONTINUE
00821    90          CONTINUE
00822             END IF
00823 *
00824             ANORM = SLANGE( '1', N, N, A, LDA, WORK )
00825             BNORM = SLANGE( '1', N, N, B, LDA, WORK )
00826 *
00827   100       CONTINUE
00828 *
00829             IF( IINFO.NE.0 ) THEN
00830                WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
00831      $            IOLDSD
00832                INFO = ABS( IINFO )
00833                RETURN
00834             END IF
00835 *
00836   110       CONTINUE
00837 *
00838 *           Call SGEQR2, SORM2R, and SGGHRD to compute H, T, U, and V
00839 *
00840             CALL SLACPY( ' ', N, N, A, LDA, H, LDA )
00841             CALL SLACPY( ' ', N, N, B, LDA, T, LDA )
00842             NTEST = 1
00843             RESULT( 1 ) = ULPINV
00844 *
00845             CALL SGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
00846             IF( IINFO.NE.0 ) THEN
00847                WRITE( NOUNIT, FMT = 9999 )'SGEQR2', IINFO, N, JTYPE,
00848      $            IOLDSD
00849                INFO = ABS( IINFO )
00850                GO TO 210
00851             END IF
00852 *
00853             CALL SORM2R( 'L', 'T', N, N, N, T, LDA, WORK, H, LDA,
00854      $                   WORK( N+1 ), IINFO )
00855             IF( IINFO.NE.0 ) THEN
00856                WRITE( NOUNIT, FMT = 9999 )'SORM2R', IINFO, N, JTYPE,
00857      $            IOLDSD
00858                INFO = ABS( IINFO )
00859                GO TO 210
00860             END IF
00861 *
00862             CALL SLASET( 'Full', N, N, ZERO, ONE, U, LDU )
00863             CALL SORM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU,
00864      $                   WORK( N+1 ), IINFO )
00865             IF( IINFO.NE.0 ) THEN
00866                WRITE( NOUNIT, FMT = 9999 )'SORM2R', IINFO, N, JTYPE,
00867      $            IOLDSD
00868                INFO = ABS( IINFO )
00869                GO TO 210
00870             END IF
00871 *
00872             CALL SGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V,
00873      $                   LDU, IINFO )
00874             IF( IINFO.NE.0 ) THEN
00875                WRITE( NOUNIT, FMT = 9999 )'SGGHRD', IINFO, N, JTYPE,
00876      $            IOLDSD
00877                INFO = ABS( IINFO )
00878                GO TO 210
00879             END IF
00880             NTEST = 4
00881 *
00882 *           Do tests 1--4
00883 *
00884             CALL SGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
00885      $                   RESULT( 1 ) )
00886             CALL SGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
00887      $                   RESULT( 2 ) )
00888             CALL SGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
00889      $                   RESULT( 3 ) )
00890             CALL SGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
00891      $                   RESULT( 4 ) )
00892 *
00893 *           Call SHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
00894 *
00895 *           Compute T1 and UZ
00896 *
00897 *           Eigenvalues only
00898 *
00899             CALL SLACPY( ' ', N, N, H, LDA, S2, LDA )
00900             CALL SLACPY( ' ', N, N, T, LDA, P2, LDA )
00901             NTEST = 5
00902             RESULT( 5 ) = ULPINV
00903 *
00904             CALL SHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
00905      $                   ALPHR3, ALPHI3, BETA3, Q, LDU, Z, LDU, WORK,
00906      $                   LWORK, IINFO )
00907             IF( IINFO.NE.0 ) THEN
00908                WRITE( NOUNIT, FMT = 9999 )'SHGEQZ(E)', IINFO, N, JTYPE,
00909      $            IOLDSD
00910                INFO = ABS( IINFO )
00911                GO TO 210
00912             END IF
00913 *
00914 *           Eigenvalues and Full Schur Form
00915 *
00916             CALL SLACPY( ' ', N, N, H, LDA, S2, LDA )
00917             CALL SLACPY( ' ', N, N, T, LDA, P2, LDA )
00918 *
00919             CALL SHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
00920      $                   ALPHR1, ALPHI1, BETA1, Q, LDU, Z, LDU, WORK,
00921      $                   LWORK, IINFO )
00922             IF( IINFO.NE.0 ) THEN
00923                WRITE( NOUNIT, FMT = 9999 )'SHGEQZ(S)', IINFO, N, JTYPE,
00924      $            IOLDSD
00925                INFO = ABS( IINFO )
00926                GO TO 210
00927             END IF
00928 *
00929 *           Eigenvalues, Schur Form, and Schur Vectors
00930 *
00931             CALL SLACPY( ' ', N, N, H, LDA, S1, LDA )
00932             CALL SLACPY( ' ', N, N, T, LDA, P1, LDA )
00933 *
00934             CALL SHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
00935      $                   ALPHR1, ALPHI1, BETA1, Q, LDU, Z, LDU, WORK,
00936      $                   LWORK, IINFO )
00937             IF( IINFO.NE.0 ) THEN
00938                WRITE( NOUNIT, FMT = 9999 )'SHGEQZ(V)', IINFO, N, JTYPE,
00939      $            IOLDSD
00940                INFO = ABS( IINFO )
00941                GO TO 210
00942             END IF
00943 *
00944             NTEST = 8
00945 *
00946 *           Do Tests 5--8
00947 *
00948             CALL SGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
00949      $                   RESULT( 5 ) )
00950             CALL SGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
00951      $                   RESULT( 6 ) )
00952             CALL SGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
00953      $                   RESULT( 7 ) )
00954             CALL SGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
00955      $                   RESULT( 8 ) )
00956 *
00957 *           Compute the Left and Right Eigenvectors of (S1,P1)
00958 *
00959 *           9: Compute the left eigenvector Matrix without
00960 *              back transforming:
00961 *
00962             NTEST = 9
00963             RESULT( 9 ) = ULPINV
00964 *
00965 *           To test "SELECT" option, compute half of the eigenvectors
00966 *           in one call, and half in another
00967 *
00968             I1 = N / 2
00969             DO 120 J = 1, I1
00970                LLWORK( J ) = .TRUE.
00971   120       CONTINUE
00972             DO 130 J = I1 + 1, N
00973                LLWORK( J ) = .FALSE.
00974   130       CONTINUE
00975 *
00976             CALL STGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
00977      $                   LDU, DUMMA, LDU, N, IN, WORK, IINFO )
00978             IF( IINFO.NE.0 ) THEN
00979                WRITE( NOUNIT, FMT = 9999 )'STGEVC(L,S1)', IINFO, N,
00980      $            JTYPE, IOLDSD
00981                INFO = ABS( IINFO )
00982                GO TO 210
00983             END IF
00984 *
00985             I1 = IN
00986             DO 140 J = 1, I1
00987                LLWORK( J ) = .FALSE.
00988   140       CONTINUE
00989             DO 150 J = I1 + 1, N
00990                LLWORK( J ) = .TRUE.
00991   150       CONTINUE
00992 *
00993             CALL STGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
00994      $                   EVECTL( 1, I1+1 ), LDU, DUMMA, LDU, N, IN,
00995      $                   WORK, IINFO )
00996             IF( IINFO.NE.0 ) THEN
00997                WRITE( NOUNIT, FMT = 9999 )'STGEVC(L,S2)', IINFO, N,
00998      $            JTYPE, IOLDSD
00999                INFO = ABS( IINFO )
01000                GO TO 210
01001             END IF
01002 *
01003             CALL SGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
01004      $                   ALPHR1, ALPHI1, BETA1, WORK, DUMMA( 1 ) )
01005             RESULT( 9 ) = DUMMA( 1 )
01006             IF( DUMMA( 2 ).GT.THRSHN ) THEN
01007                WRITE( NOUNIT, FMT = 9998 )'Left', 'STGEVC(HOWMNY=S)',
01008      $            DUMMA( 2 ), N, JTYPE, IOLDSD
01009             END IF
01010 *
01011 *           10: Compute the left eigenvector Matrix with
01012 *               back transforming:
01013 *
01014             NTEST = 10
01015             RESULT( 10 ) = ULPINV
01016             CALL SLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
01017             CALL STGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
01018      $                   LDU, DUMMA, LDU, N, IN, WORK, IINFO )
01019             IF( IINFO.NE.0 ) THEN
01020                WRITE( NOUNIT, FMT = 9999 )'STGEVC(L,B)', IINFO, N,
01021      $            JTYPE, IOLDSD
01022                INFO = ABS( IINFO )
01023                GO TO 210
01024             END IF
01025 *
01026             CALL SGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHR1,
01027      $                   ALPHI1, BETA1, WORK, DUMMA( 1 ) )
01028             RESULT( 10 ) = DUMMA( 1 )
01029             IF( DUMMA( 2 ).GT.THRSHN ) THEN
01030                WRITE( NOUNIT, FMT = 9998 )'Left', 'STGEVC(HOWMNY=B)',
01031      $            DUMMA( 2 ), N, JTYPE, IOLDSD
01032             END IF
01033 *
01034 *           11: Compute the right eigenvector Matrix without
01035 *               back transforming:
01036 *
01037             NTEST = 11
01038             RESULT( 11 ) = ULPINV
01039 *
01040 *           To test "SELECT" option, compute half of the eigenvectors
01041 *           in one call, and half in another
01042 *
01043             I1 = N / 2
01044             DO 160 J = 1, I1
01045                LLWORK( J ) = .TRUE.
01046   160       CONTINUE
01047             DO 170 J = I1 + 1, N
01048                LLWORK( J ) = .FALSE.
01049   170       CONTINUE
01050 *
01051             CALL STGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, DUMMA,
01052      $                   LDU, EVECTR, LDU, N, IN, WORK, IINFO )
01053             IF( IINFO.NE.0 ) THEN
01054                WRITE( NOUNIT, FMT = 9999 )'STGEVC(R,S1)', IINFO, N,
01055      $            JTYPE, IOLDSD
01056                INFO = ABS( IINFO )
01057                GO TO 210
01058             END IF
01059 *
01060             I1 = IN
01061             DO 180 J = 1, I1
01062                LLWORK( J ) = .FALSE.
01063   180       CONTINUE
01064             DO 190 J = I1 + 1, N
01065                LLWORK( J ) = .TRUE.
01066   190       CONTINUE
01067 *
01068             CALL STGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, DUMMA,
01069      $                   LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
01070      $                   IINFO )
01071             IF( IINFO.NE.0 ) THEN
01072                WRITE( NOUNIT, FMT = 9999 )'STGEVC(R,S2)', IINFO, N,
01073      $            JTYPE, IOLDSD
01074                INFO = ABS( IINFO )
01075                GO TO 210
01076             END IF
01077 *
01078             CALL SGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
01079      $                   ALPHR1, ALPHI1, BETA1, WORK, DUMMA( 1 ) )
01080             RESULT( 11 ) = DUMMA( 1 )
01081             IF( DUMMA( 2 ).GT.THRESH ) THEN
01082                WRITE( NOUNIT, FMT = 9998 )'Right', 'STGEVC(HOWMNY=S)',
01083      $            DUMMA( 2 ), N, JTYPE, IOLDSD
01084             END IF
01085 *
01086 *           12: Compute the right eigenvector Matrix with
01087 *               back transforming:
01088 *
01089             NTEST = 12
01090             RESULT( 12 ) = ULPINV
01091             CALL SLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
01092             CALL STGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, DUMMA,
01093      $                   LDU, EVECTR, LDU, N, IN, WORK, IINFO )
01094             IF( IINFO.NE.0 ) THEN
01095                WRITE( NOUNIT, FMT = 9999 )'STGEVC(R,B)', IINFO, N,
01096      $            JTYPE, IOLDSD
01097                INFO = ABS( IINFO )
01098                GO TO 210
01099             END IF
01100 *
01101             CALL SGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
01102      $                   ALPHR1, ALPHI1, BETA1, WORK, DUMMA( 1 ) )
01103             RESULT( 12 ) = DUMMA( 1 )
01104             IF( DUMMA( 2 ).GT.THRESH ) THEN
01105                WRITE( NOUNIT, FMT = 9998 )'Right', 'STGEVC(HOWMNY=B)',
01106      $            DUMMA( 2 ), N, JTYPE, IOLDSD
01107             END IF
01108 *
01109 *           Tests 13--15 are done only on request
01110 *
01111             IF( TSTDIF ) THEN
01112 *
01113 *              Do Tests 13--14
01114 *
01115                CALL SGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
01116      $                      WORK, RESULT( 13 ) )
01117                CALL SGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
01118      $                      WORK, RESULT( 14 ) )
01119 *
01120 *              Do Test 15
01121 *
01122                TEMP1 = ZERO
01123                TEMP2 = ZERO
01124                DO 200 J = 1, N
01125                   TEMP1 = MAX( TEMP1, ABS( ALPHR1( J )-ALPHR3( J ) )+
01126      $                    ABS( ALPHI1( J )-ALPHI3( J ) ) )
01127                   TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
01128   200          CONTINUE
01129 *
01130                TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
01131                TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
01132                RESULT( 15 ) = MAX( TEMP1, TEMP2 )
01133                NTEST = 15
01134             ELSE
01135                RESULT( 13 ) = ZERO
01136                RESULT( 14 ) = ZERO
01137                RESULT( 15 ) = ZERO
01138                NTEST = 12
01139             END IF
01140 *
01141 *           End of Loop -- Check for RESULT(j) > THRESH
01142 *
01143   210       CONTINUE
01144 *
01145             NTESTT = NTESTT + NTEST
01146 *
01147 *           Print out tests which fail.
01148 *
01149             DO 220 JR = 1, NTEST
01150                IF( RESULT( JR ).GE.THRESH ) THEN
01151 *
01152 *                 If this is the first test to fail,
01153 *                 print a header to the data file.
01154 *
01155                   IF( NERRS.EQ.0 ) THEN
01156                      WRITE( NOUNIT, FMT = 9997 )'SGG'
01157 *
01158 *                    Matrix types
01159 *
01160                      WRITE( NOUNIT, FMT = 9996 )
01161                      WRITE( NOUNIT, FMT = 9995 )
01162                      WRITE( NOUNIT, FMT = 9994 )'Orthogonal'
01163 *
01164 *                    Tests performed
01165 *
01166                      WRITE( NOUNIT, FMT = 9993 )'orthogonal', '''',
01167      $                  'transpose', ( '''', J = 1, 10 )
01168 *
01169                   END IF
01170                   NERRS = NERRS + 1
01171                   IF( RESULT( JR ).LT.10000.0 ) THEN
01172                      WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
01173      $                  RESULT( JR )
01174                   ELSE
01175                      WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
01176      $                  RESULT( JR )
01177                   END IF
01178                END IF
01179   220       CONTINUE
01180 *
01181   230    CONTINUE
01182   240 CONTINUE
01183 *
01184 *     Summary
01185 *
01186       CALL SLASUM( 'SGG', NOUNIT, NERRS, NTESTT )
01187       RETURN
01188 *
01189  9999 FORMAT( ' SCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
01190      $      I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
01191 *
01192  9998 FORMAT( ' SCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
01193      $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
01194      $      'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
01195      $      ')' )
01196 *
01197  9997 FORMAT( / 1X, A3, ' -- Real Generalized eigenvalue problem' )
01198 *
01199  9996 FORMAT( ' Matrix types (see SCHKGG for details): ' )
01200 *
01201  9995 FORMAT( ' Special Matrices:', 23X,
01202      $      '(J''=transposed Jordan block)',
01203      $      / '   1=(0,0)  2=(I,0)  3=(0,I)  4=(I,I)  5=(J'',J'')  ',
01204      $      '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices:  ( ',
01205      $      'D=diag(0,1,2,...) )', / '   7=(D,I)   9=(large*D, small*I',
01206      $      ')  11=(large*I, small*D)  13=(large*D, large*I)', /
01207      $      '   8=(I,D)  10=(small*D, large*I)  12=(small*I, large*D) ',
01208      $      ' 14=(small*D, small*I)', / '  15=(D, reversed D)' )
01209  9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
01210      $      / '  16=Transposed Jordan Blocks             19=geometric ',
01211      $      'alpha, beta=0,1', / '  17=arithm. alpha&beta             ',
01212      $      '      20=arithmetic alpha, beta=0,1', / '  18=clustered ',
01213      $      'alpha, beta=0,1            21=random alpha, beta=0,1',
01214      $      / ' Large & Small Matrices:', / '  22=(large, small)   ',
01215      $      '23=(small,large)    24=(small,small)    25=(large,large)',
01216      $      / '  26=random O(1) matrices.' )
01217 *
01218  9993 FORMAT( / ' Tests performed:   (H is Hessenberg, S is Schur, B, ',
01219      $      'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
01220      $      ', l and r are the', / 20X,
01221      $      'appropriate left and right eigenvectors, resp., a is',
01222      $      / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
01223      $      / ' 1 = | A - U H V', A,
01224      $      ' | / ( |A| n ulp )      2 = | B - U T V', A,
01225      $      ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
01226      $      ' | / ( n ulp )             4 = | I - VV', A,
01227      $      ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
01228      $      ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
01229      $      ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
01230      $      ' | / ( n ulp )             8 = | I - ZZ', A,
01231      $      ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
01232      $      ' l | / const.  10 = max | ( b H - a T )', A,
01233      $      ' l | / const.', /
01234      $      ' 11= max | ( b S - a P ) r | / const.   12 = max | ( b H',
01235      $      ' - a T ) r | / const.', / 1X )
01236 *
01237  9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
01238      $      4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
01239  9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
01240      $      4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
01241 *
01242 *     End of SCHKGG
01243 *
01244       END
 All Files Functions