LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
schkhs.f
Go to the documentation of this file.
1 *> \brief \b SCHKHS
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE SCHKHS( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, H, T1, T2, U, LDU, Z, UZ, WR1,
13 * WI1, WR3, WI3, EVECTL, EVECTR, EVECTY, EVECTX,
14 * UU, TAU, WORK, NWORK, IWORK, SELECT, RESULT,
15 * INFO )
16 *
17 * .. Scalar Arguments ..
18 * INTEGER INFO, LDA, LDU, NOUNIT, NSIZES, NTYPES, NWORK
19 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * ), SELECT( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24 * REAL A( LDA, * ), EVECTL( LDU, * ),
25 * $ EVECTR( LDU, * ), EVECTX( LDU, * ),
26 * $ EVECTY( LDU, * ), H( LDA, * ), RESULT( 14 ),
27 * $ T1( LDA, * ), T2( LDA, * ), TAU( * ),
28 * $ U( LDU, * ), UU( LDU, * ), UZ( LDU, * ),
29 * $ WI1( * ), WI3( * ), WORK( * ), WR1( * ),
30 * $ WR3( * ), Z( LDU, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SCHKHS checks the nonsymmetric eigenvalue problem routines.
40 *>
41 *> SGEHRD factors A as U H U' , where ' means transpose,
42 *> H is hessenberg, and U is an orthogonal matrix.
43 *>
44 *> SORGHR generates the orthogonal matrix U.
45 *>
46 *> SORMHR multiplies a matrix by the orthogonal matrix U.
47 *>
48 *> SHSEQR factors H as Z T Z' , where Z is orthogonal and
49 *> T is "quasi-triangular", and the eigenvalue vector W.
50 *>
51 *> STREVC computes the left and right eigenvector matrices
52 *> L and R for T.
53 *>
54 *> SHSEIN computes the left and right eigenvector matrices
55 *> Y and X for H, using inverse iteration.
56 *>
57 *> When SCHKHS is called, a number of matrix "sizes" ("n's") and a
58 *> number of matrix "types" are specified. For each size ("n")
59 *> and each type of matrix, one matrix will be generated and used
60 *> to test the nonsymmetric eigenroutines. For each matrix, 14
61 *> tests will be performed:
62 *>
63 *> (1) | A - U H U**T | / ( |A| n ulp )
64 *>
65 *> (2) | I - UU**T | / ( n ulp )
66 *>
67 *> (3) | H - Z T Z**T | / ( |H| n ulp )
68 *>
69 *> (4) | I - ZZ**T | / ( n ulp )
70 *>
71 *> (5) | A - UZ H (UZ)**T | / ( |A| n ulp )
72 *>
73 *> (6) | I - UZ (UZ)**T | / ( n ulp )
74 *>
75 *> (7) | T(Z computed) - T(Z not computed) | / ( |T| ulp )
76 *>
77 *> (8) | W(Z computed) - W(Z not computed) | / ( |W| ulp )
78 *>
79 *> (9) | TR - RW | / ( |T| |R| ulp )
80 *>
81 *> (10) | L**H T - W**H L | / ( |T| |L| ulp )
82 *>
83 *> (11) | HX - XW | / ( |H| |X| ulp )
84 *>
85 *> (12) | Y**H H - W**H Y | / ( |H| |Y| ulp )
86 *>
87 *> (13) | AX - XW | / ( |A| |X| ulp )
88 *>
89 *> (14) | Y**H A - W**H Y | / ( |A| |Y| ulp )
90 *>
91 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
92 *> each element NN(j) specifies one size.
93 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
94 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
95 *> Currently, the list of possible types is:
96 *>
97 *> (1) The zero matrix.
98 *> (2) The identity matrix.
99 *> (3) A (transposed) Jordan block, with 1's on the diagonal.
100 *>
101 *> (4) A diagonal matrix with evenly spaced entries
102 *> 1, ..., ULP and random signs.
103 *> (ULP = (first number larger than 1) - 1 )
104 *> (5) A diagonal matrix with geometrically spaced entries
105 *> 1, ..., ULP and random signs.
106 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
107 *> and random signs.
108 *>
109 *> (7) Same as (4), but multiplied by SQRT( overflow threshold )
110 *> (8) Same as (4), but multiplied by SQRT( underflow threshold )
111 *>
112 *> (9) A matrix of the form U' T U, where U is orthogonal and
113 *> T has evenly spaced entries 1, ..., ULP with random signs
114 *> on the diagonal and random O(1) entries in the upper
115 *> triangle.
116 *>
117 *> (10) A matrix of the form U' T U, where U is orthogonal and
118 *> T has geometrically spaced entries 1, ..., ULP with random
119 *> signs on the diagonal and random O(1) entries in the upper
120 *> triangle.
121 *>
122 *> (11) A matrix of the form U' T U, where U is orthogonal and
123 *> T has "clustered" entries 1, ULP,..., ULP with random
124 *> signs on the diagonal and random O(1) entries in the upper
125 *> triangle.
126 *>
127 *> (12) A matrix of the form U' T U, where U is orthogonal and
128 *> T has real or complex conjugate paired eigenvalues randomly
129 *> chosen from ( ULP, 1 ) and random O(1) entries in the upper
130 *> triangle.
131 *>
132 *> (13) A matrix of the form X' T X, where X has condition
133 *> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
134 *> with random signs on the diagonal and random O(1) entries
135 *> in the upper triangle.
136 *>
137 *> (14) A matrix of the form X' T X, where X has condition
138 *> SQRT( ULP ) and T has geometrically spaced entries
139 *> 1, ..., ULP with random signs on the diagonal and random
140 *> O(1) entries in the upper triangle.
141 *>
142 *> (15) A matrix of the form X' T X, where X has condition
143 *> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
144 *> with random signs on the diagonal and random O(1) entries
145 *> in the upper triangle.
146 *>
147 *> (16) A matrix of the form X' T X, where X has condition
148 *> SQRT( ULP ) and T has real or complex conjugate paired
149 *> eigenvalues randomly chosen from ( ULP, 1 ) and random
150 *> O(1) entries in the upper triangle.
151 *>
152 *> (17) Same as (16), but multiplied by SQRT( overflow threshold )
153 *> (18) Same as (16), but multiplied by SQRT( underflow threshold )
154 *>
155 *> (19) Nonsymmetric matrix with random entries chosen from (-1,1).
156 *> (20) Same as (19), but multiplied by SQRT( overflow threshold )
157 *> (21) Same as (19), but multiplied by SQRT( underflow threshold )
158 *> \endverbatim
159 *
160 * Arguments:
161 * ==========
162 *
163 *> \verbatim
164 *> NSIZES - INTEGER
165 *> The number of sizes of matrices to use. If it is zero,
166 *> SCHKHS does nothing. It must be at least zero.
167 *> Not modified.
168 *>
169 *> NN - INTEGER array, dimension (NSIZES)
170 *> An array containing the sizes to be used for the matrices.
171 *> Zero values will be skipped. The values must be at least
172 *> zero.
173 *> Not modified.
174 *>
175 *> NTYPES - INTEGER
176 *> The number of elements in DOTYPE. If it is zero, SCHKHS
177 *> does nothing. It must be at least zero. If it is MAXTYP+1
178 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
179 *> defined, which is to use whatever matrix is in A. This
180 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
181 *> DOTYPE(MAXTYP+1) is .TRUE. .
182 *> Not modified.
183 *>
184 *> DOTYPE - LOGICAL array, dimension (NTYPES)
185 *> If DOTYPE(j) is .TRUE., then for each size in NN a
186 *> matrix of that size and of type j will be generated.
187 *> If NTYPES is smaller than the maximum number of types
188 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
189 *> MAXTYP will not be generated. If NTYPES is larger
190 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
191 *> will be ignored.
192 *> Not modified.
193 *>
194 *> ISEED - INTEGER array, dimension (4)
195 *> On entry ISEED specifies the seed of the random number
196 *> generator. The array elements should be between 0 and 4095;
197 *> if not they will be reduced mod 4096. Also, ISEED(4) must
198 *> be odd. The random number generator uses a linear
199 *> congruential sequence limited to small integers, and so
200 *> should produce machine independent random numbers. The
201 *> values of ISEED are changed on exit, and can be used in the
202 *> next call to SCHKHS to continue the same random number
203 *> sequence.
204 *> Modified.
205 *>
206 *> THRESH - REAL
207 *> A test will count as "failed" if the "error", computed as
208 *> described above, exceeds THRESH. Note that the error
209 *> is scaled to be O(1), so THRESH should be a reasonably
210 *> small multiple of 1, e.g., 10 or 100. In particular,
211 *> it should not depend on the precision (single vs. double)
212 *> or the size of the matrix. It must be at least zero.
213 *> Not modified.
214 *>
215 *> NOUNIT - INTEGER
216 *> The FORTRAN unit number for printing out error messages
217 *> (e.g., if a routine returns IINFO not equal to 0.)
218 *> Not modified.
219 *>
220 *> A - REAL array, dimension (LDA,max(NN))
221 *> Used to hold the matrix whose eigenvalues are to be
222 *> computed. On exit, A contains the last matrix actually
223 *> used.
224 *> Modified.
225 *>
226 *> LDA - INTEGER
227 *> The leading dimension of A, H, T1 and T2. It must be at
228 *> least 1 and at least max( NN ).
229 *> Not modified.
230 *>
231 *> H - REAL array, dimension (LDA,max(NN))
232 *> The upper hessenberg matrix computed by SGEHRD. On exit,
233 *> H contains the Hessenberg form of the matrix in A.
234 *> Modified.
235 *>
236 *> T1 - REAL array, dimension (LDA,max(NN))
237 *> The Schur (="quasi-triangular") matrix computed by SHSEQR
238 *> if Z is computed. On exit, T1 contains the Schur form of
239 *> the matrix in A.
240 *> Modified.
241 *>
242 *> T2 - REAL array, dimension (LDA,max(NN))
243 *> The Schur matrix computed by SHSEQR when Z is not computed.
244 *> This should be identical to T1.
245 *> Modified.
246 *>
247 *> LDU - INTEGER
248 *> The leading dimension of U, Z, UZ and UU. It must be at
249 *> least 1 and at least max( NN ).
250 *> Not modified.
251 *>
252 *> U - REAL array, dimension (LDU,max(NN))
253 *> The orthogonal matrix computed by SGEHRD.
254 *> Modified.
255 *>
256 *> Z - REAL array, dimension (LDU,max(NN))
257 *> The orthogonal matrix computed by SHSEQR.
258 *> Modified.
259 *>
260 *> UZ - REAL array, dimension (LDU,max(NN))
261 *> The product of U times Z.
262 *> Modified.
263 *>
264 *> WR1 - REAL array, dimension (max(NN))
265 *> WI1 - REAL array, dimension (max(NN))
266 *> The real and imaginary parts of the eigenvalues of A,
267 *> as computed when Z is computed.
268 *> On exit, WR1 + WI1*i are the eigenvalues of the matrix in A.
269 *> Modified.
270 *>
271 *> WR3 - REAL array, dimension (max(NN))
272 *> WI3 - REAL array, dimension (max(NN))
273 *> Like WR1, WI1, these arrays contain the eigenvalues of A,
274 *> but those computed when SHSEQR only computes the
275 *> eigenvalues, i.e., not the Schur vectors and no more of the
276 *> Schur form than is necessary for computing the
277 *> eigenvalues.
278 *> Modified.
279 *>
280 *> EVECTL - REAL array, dimension (LDU,max(NN))
281 *> The (upper triangular) left eigenvector matrix for the
282 *> matrix in T1. For complex conjugate pairs, the real part
283 *> is stored in one row and the imaginary part in the next.
284 *> Modified.
285 *>
286 *> EVECTR - REAL array, dimension (LDU,max(NN))
287 *> The (upper triangular) right eigenvector matrix for the
288 *> matrix in T1. For complex conjugate pairs, the real part
289 *> is stored in one column and the imaginary part in the next.
290 *> Modified.
291 *>
292 *> EVECTY - REAL array, dimension (LDU,max(NN))
293 *> The left eigenvector matrix for the
294 *> matrix in H. For complex conjugate pairs, the real part
295 *> is stored in one row and the imaginary part in the next.
296 *> Modified.
297 *>
298 *> EVECTX - REAL array, dimension (LDU,max(NN))
299 *> The right eigenvector matrix for the
300 *> matrix in H. For complex conjugate pairs, the real part
301 *> is stored in one column and the imaginary part in the next.
302 *> Modified.
303 *>
304 *> UU - REAL array, dimension (LDU,max(NN))
305 *> Details of the orthogonal matrix computed by SGEHRD.
306 *> Modified.
307 *>
308 *> TAU - REAL array, dimension(max(NN))
309 *> Further details of the orthogonal matrix computed by SGEHRD.
310 *> Modified.
311 *>
312 *> WORK - REAL array, dimension (NWORK)
313 *> Workspace.
314 *> Modified.
315 *>
316 *> NWORK - INTEGER
317 *> The number of entries in WORK. NWORK >= 4*NN(j)*NN(j) + 2.
318 *>
319 *> IWORK - INTEGER array, dimension (max(NN))
320 *> Workspace.
321 *> Modified.
322 *>
323 *> SELECT - LOGICAL array, dimension (max(NN))
324 *> Workspace.
325 *> Modified.
326 *>
327 *> RESULT - REAL array, dimension (14)
328 *> The values computed by the fourteen tests described above.
329 *> The values are currently limited to 1/ulp, to avoid
330 *> overflow.
331 *> Modified.
332 *>
333 *> INFO - INTEGER
334 *> If 0, then everything ran OK.
335 *> -1: NSIZES < 0
336 *> -2: Some NN(j) < 0
337 *> -3: NTYPES < 0
338 *> -6: THRESH < 0
339 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
340 *> -14: LDU < 1 or LDU < NMAX.
341 *> -28: NWORK too small.
342 *> If SLATMR, SLATMS, or SLATME returns an error code, the
343 *> absolute value of it is returned.
344 *> If 1, then SHSEQR could not find all the shifts.
345 *> If 2, then the EISPACK code (for small blocks) failed.
346 *> If >2, then 30*N iterations were not enough to find an
347 *> eigenvalue or to decompose the problem.
348 *> Modified.
349 *>
350 *>-----------------------------------------------------------------------
351 *>
352 *> Some Local Variables and Parameters:
353 *> ---- ----- --------- --- ----------
354 *>
355 *> ZERO, ONE Real 0 and 1.
356 *> MAXTYP The number of types defined.
357 *> MTEST The number of tests defined: care must be taken
358 *> that (1) the size of RESULT, (2) the number of
359 *> tests actually performed, and (3) MTEST agree.
360 *> NTEST The number of tests performed on this matrix
361 *> so far. This should be less than MTEST, and
362 *> equal to it by the last test. It will be less
363 *> if any of the routines being tested indicates
364 *> that it could not compute the matrices that
365 *> would be tested.
366 *> NMAX Largest value in NN.
367 *> NMATS The number of matrices generated so far.
368 *> NERRS The number of tests which have exceeded THRESH
369 *> so far (computed by SLAFTS).
370 *> COND, CONDS,
371 *> IMODE Values to be passed to the matrix generators.
372 *> ANORM Norm of A; passed to matrix generators.
373 *>
374 *> OVFL, UNFL Overflow and underflow thresholds.
375 *> ULP, ULPINV Finest relative precision and its inverse.
376 *> RTOVFL, RTUNFL,
377 *> RTULP, RTULPI Square roots of the previous 4 values.
378 *>
379 *> The following four arrays decode JTYPE:
380 *> KTYPE(j) The general type (1-10) for type "j".
381 *> KMODE(j) The MODE value to be passed to the matrix
382 *> generator for type "j".
383 *> KMAGN(j) The order of magnitude ( O(1),
384 *> O(overflow^(1/2) ), O(underflow^(1/2) )
385 *> KCONDS(j) Selects whether CONDS is to be 1 or
386 *> 1/sqrt(ulp). (0 means irrelevant.)
387 *> \endverbatim
388 *
389 * Authors:
390 * ========
391 *
392 *> \author Univ. of Tennessee
393 *> \author Univ. of California Berkeley
394 *> \author Univ. of Colorado Denver
395 *> \author NAG Ltd.
396 *
397 *> \date November 2011
398 *
399 *> \ingroup single_eig
400 *
401 * =====================================================================
402  SUBROUTINE schkhs( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
403  $ nounit, a, lda, h, t1, t2, u, ldu, z, uz, wr1,
404  $ wi1, wr3, wi3, evectl, evectr, evecty, evectx,
405  $ uu, tau, work, nwork, iwork, SELECT, result,
406  $ info )
407 *
408 * -- LAPACK test routine (version 3.4.0) --
409 * -- LAPACK is a software package provided by Univ. of Tennessee, --
410 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
411 * November 2011
412 *
413 * .. Scalar Arguments ..
414  INTEGER info, lda, ldu, nounit, nsizes, ntypes, nwork
415  REAL thresh
416 * ..
417 * .. Array Arguments ..
418  LOGICAL dotype( * ), select( * )
419  INTEGER iseed( 4 ), iwork( * ), nn( * )
420  REAL a( lda, * ), evectl( ldu, * ),
421  $ evectr( ldu, * ), evectx( ldu, * ),
422  $ evecty( ldu, * ), h( lda, * ), result( 14 ),
423  $ t1( lda, * ), t2( lda, * ), tau( * ),
424  $ u( ldu, * ), uu( ldu, * ), uz( ldu, * ),
425  $ wi1( * ), wi3( * ), work( * ), wr1( * ),
426  $ wr3( * ), z( ldu, * )
427 * ..
428 *
429 * =====================================================================
430 *
431 * .. Parameters ..
432  REAL zero, one
433  parameter( zero = 0.0, one = 1.0 )
434  INTEGER maxtyp
435  parameter( maxtyp = 21 )
436 * ..
437 * .. Local Scalars ..
438  LOGICAL badnn, match
439  INTEGER i, ihi, iinfo, ilo, imode, in, itype, j, jcol,
440  $ jj, jsize, jtype, k, mtypes, n, n1, nerrs,
441  $ nmats, nmax, nselc, nselr, ntest, ntestt
442  REAL aninv, anorm, cond, conds, ovfl, rtovfl, rtulp,
443  $ rtulpi, rtunfl, temp1, temp2, ulp, ulpinv, unfl
444 * ..
445 * .. Local Arrays ..
446  CHARACTER adumma( 1 )
447  INTEGER idumma( 1 ), ioldsd( 4 ), kconds( maxtyp ),
448  $ kmagn( maxtyp ), kmode( maxtyp ),
449  $ ktype( maxtyp )
450  REAL dumma( 6 )
451 * ..
452 * .. External Functions ..
453  REAL slamch
454  EXTERNAL slamch
455 * ..
456 * .. External Subroutines ..
457  EXTERNAL scopy, sgehrd, sgemm, sget10, sget22, shsein,
458  $ shseqr, shst01, slabad, slacpy, slafts, slaset,
459  $ slasum, slatme, slatmr, slatms, sorghr, sormhr,
460  $ strevc, xerbla
461 * ..
462 * .. Intrinsic Functions ..
463  INTRINSIC abs, max, min, REAL, sqrt
464 * ..
465 * .. Data statements ..
466  DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
467  DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
468  $ 3, 1, 2, 3 /
469  DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
470  $ 1, 5, 5, 5, 4, 3, 1 /
471  DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
472 * ..
473 * .. Executable Statements ..
474 *
475 * Check for errors
476 *
477  ntestt = 0
478  info = 0
479 *
480  badnn = .false.
481  nmax = 0
482  DO 10 j = 1, nsizes
483  nmax = max( nmax, nn( j ) )
484  IF( nn( j ).LT.0 )
485  $ badnn = .true.
486  10 continue
487 *
488 * Check for errors
489 *
490  IF( nsizes.LT.0 ) THEN
491  info = -1
492  ELSE IF( badnn ) THEN
493  info = -2
494  ELSE IF( ntypes.LT.0 ) THEN
495  info = -3
496  ELSE IF( thresh.LT.zero ) THEN
497  info = -6
498  ELSE IF( lda.LE.1 .OR. lda.LT.nmax ) THEN
499  info = -9
500  ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax ) THEN
501  info = -14
502  ELSE IF( 4*nmax*nmax+2.GT.nwork ) THEN
503  info = -28
504  END IF
505 *
506  IF( info.NE.0 ) THEN
507  CALL xerbla( 'SCHKHS', -info )
508  return
509  END IF
510 *
511 * Quick return if possible
512 *
513  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
514  $ return
515 *
516 * More important constants
517 *
518  unfl = slamch( 'Safe minimum' )
519  ovfl = slamch( 'Overflow' )
520  CALL slabad( unfl, ovfl )
521  ulp = slamch( 'Epsilon' )*slamch( 'Base' )
522  ulpinv = one / ulp
523  rtunfl = sqrt( unfl )
524  rtovfl = sqrt( ovfl )
525  rtulp = sqrt( ulp )
526  rtulpi = one / rtulp
527 *
528 * Loop over sizes, types
529 *
530  nerrs = 0
531  nmats = 0
532 *
533  DO 270 jsize = 1, nsizes
534  n = nn( jsize )
535  IF( n.EQ.0 )
536  $ go to 270
537  n1 = max( 1, n )
538  aninv = one / REAL( n1 )
539 *
540  IF( nsizes.NE.1 ) THEN
541  mtypes = min( maxtyp, ntypes )
542  ELSE
543  mtypes = min( maxtyp+1, ntypes )
544  END IF
545 *
546  DO 260 jtype = 1, mtypes
547  IF( .NOT.dotype( jtype ) )
548  $ go to 260
549  nmats = nmats + 1
550  ntest = 0
551 *
552 * Save ISEED in case of an error.
553 *
554  DO 20 j = 1, 4
555  ioldsd( j ) = iseed( j )
556  20 continue
557 *
558 * Initialize RESULT
559 *
560  DO 30 j = 1, 14
561  result( j ) = zero
562  30 continue
563 *
564 * Compute "A"
565 *
566 * Control parameters:
567 *
568 * KMAGN KCONDS KMODE KTYPE
569 * =1 O(1) 1 clustered 1 zero
570 * =2 large large clustered 2 identity
571 * =3 small exponential Jordan
572 * =4 arithmetic diagonal, (w/ eigenvalues)
573 * =5 random log symmetric, w/ eigenvalues
574 * =6 random general, w/ eigenvalues
575 * =7 random diagonal
576 * =8 random symmetric
577 * =9 random general
578 * =10 random triangular
579 *
580  IF( mtypes.GT.maxtyp )
581  $ go to 100
582 *
583  itype = ktype( jtype )
584  imode = kmode( jtype )
585 *
586 * Compute norm
587 *
588  go to( 40, 50, 60 )kmagn( jtype )
589 *
590  40 continue
591  anorm = one
592  go to 70
593 *
594  50 continue
595  anorm = ( rtovfl*ulp )*aninv
596  go to 70
597 *
598  60 continue
599  anorm = rtunfl*n*ulpinv
600  go to 70
601 *
602  70 continue
603 *
604  CALL slaset( 'Full', lda, n, zero, zero, a, lda )
605  iinfo = 0
606  cond = ulpinv
607 *
608 * Special Matrices
609 *
610  IF( itype.EQ.1 ) THEN
611 *
612 * Zero
613 *
614  iinfo = 0
615 *
616  ELSE IF( itype.EQ.2 ) THEN
617 *
618 * Identity
619 *
620  DO 80 jcol = 1, n
621  a( jcol, jcol ) = anorm
622  80 continue
623 *
624  ELSE IF( itype.EQ.3 ) THEN
625 *
626 * Jordan Block
627 *
628  DO 90 jcol = 1, n
629  a( jcol, jcol ) = anorm
630  IF( jcol.GT.1 )
631  $ a( jcol, jcol-1 ) = one
632  90 continue
633 *
634  ELSE IF( itype.EQ.4 ) THEN
635 *
636 * Diagonal Matrix, [Eigen]values Specified
637 *
638  CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
639  $ anorm, 0, 0, 'N', a, lda, work( n+1 ),
640  $ iinfo )
641 *
642  ELSE IF( itype.EQ.5 ) THEN
643 *
644 * Symmetric, eigenvalues specified
645 *
646  CALL slatms( n, n, 'S', iseed, 'S', work, imode, cond,
647  $ anorm, n, n, 'N', a, lda, work( n+1 ),
648  $ iinfo )
649 *
650  ELSE IF( itype.EQ.6 ) THEN
651 *
652 * General, eigenvalues specified
653 *
654  IF( kconds( jtype ).EQ.1 ) THEN
655  conds = one
656  ELSE IF( kconds( jtype ).EQ.2 ) THEN
657  conds = rtulpi
658  ELSE
659  conds = zero
660  END IF
661 *
662  adumma( 1 ) = ' '
663  CALL slatme( n, 'S', iseed, work, imode, cond, one,
664  $ adumma, 'T', 'T', 'T', work( n+1 ), 4,
665  $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
666  $ iinfo )
667 *
668  ELSE IF( itype.EQ.7 ) THEN
669 *
670 * Diagonal, random eigenvalues
671 *
672  CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
673  $ 'T', 'N', work( n+1 ), 1, one,
674  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
675  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
676 *
677  ELSE IF( itype.EQ.8 ) THEN
678 *
679 * Symmetric, random eigenvalues
680 *
681  CALL slatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
682  $ 'T', 'N', work( n+1 ), 1, one,
683  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
684  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
685 *
686  ELSE IF( itype.EQ.9 ) THEN
687 *
688 * General, random eigenvalues
689 *
690  CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
691  $ 'T', 'N', work( n+1 ), 1, one,
692  $ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
693  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
694 *
695  ELSE IF( itype.EQ.10 ) THEN
696 *
697 * Triangular, random eigenvalues
698 *
699  CALL slatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
700  $ 'T', 'N', work( n+1 ), 1, one,
701  $ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
702  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
703 *
704  ELSE
705 *
706  iinfo = 1
707  END IF
708 *
709  IF( iinfo.NE.0 ) THEN
710  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n, jtype,
711  $ ioldsd
712  info = abs( iinfo )
713  return
714  END IF
715 *
716  100 continue
717 *
718 * Call SGEHRD to compute H and U, do tests.
719 *
720  CALL slacpy( ' ', n, n, a, lda, h, lda )
721 *
722  ntest = 1
723 *
724  ilo = 1
725  ihi = n
726 *
727  CALL sgehrd( n, ilo, ihi, h, lda, work, work( n+1 ),
728  $ nwork-n, iinfo )
729 *
730  IF( iinfo.NE.0 ) THEN
731  result( 1 ) = ulpinv
732  WRITE( nounit, fmt = 9999 )'SGEHRD', iinfo, n, jtype,
733  $ ioldsd
734  info = abs( iinfo )
735  go to 250
736  END IF
737 *
738  DO 120 j = 1, n - 1
739  uu( j+1, j ) = zero
740  DO 110 i = j + 2, n
741  u( i, j ) = h( i, j )
742  uu( i, j ) = h( i, j )
743  h( i, j ) = zero
744  110 continue
745  120 continue
746  CALL scopy( n-1, work, 1, tau, 1 )
747  CALL sorghr( n, ilo, ihi, u, ldu, work, work( n+1 ),
748  $ nwork-n, iinfo )
749  ntest = 2
750 *
751  CALL shst01( n, ilo, ihi, a, lda, h, lda, u, ldu, work,
752  $ nwork, result( 1 ) )
753 *
754 * Call SHSEQR to compute T1, T2 and Z, do tests.
755 *
756 * Eigenvalues only (WR3,WI3)
757 *
758  CALL slacpy( ' ', n, n, h, lda, t2, lda )
759  ntest = 3
760  result( 3 ) = ulpinv
761 *
762  CALL shseqr( 'E', 'N', n, ilo, ihi, t2, lda, wr3, wi3, uz,
763  $ ldu, work, nwork, iinfo )
764  IF( iinfo.NE.0 ) THEN
765  WRITE( nounit, fmt = 9999 )'SHSEQR(E)', iinfo, n, jtype,
766  $ ioldsd
767  IF( iinfo.LE.n+2 ) THEN
768  info = abs( iinfo )
769  go to 250
770  END IF
771  END IF
772 *
773 * Eigenvalues (WR1,WI1) and Full Schur Form (T2)
774 *
775  CALL slacpy( ' ', n, n, h, lda, t2, lda )
776 *
777  CALL shseqr( 'S', 'N', n, ilo, ihi, t2, lda, wr1, wi1, uz,
778  $ ldu, work, nwork, iinfo )
779  IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
780  WRITE( nounit, fmt = 9999 )'SHSEQR(S)', iinfo, n, jtype,
781  $ ioldsd
782  info = abs( iinfo )
783  go to 250
784  END IF
785 *
786 * Eigenvalues (WR1,WI1), Schur Form (T1), and Schur vectors
787 * (UZ)
788 *
789  CALL slacpy( ' ', n, n, h, lda, t1, lda )
790  CALL slacpy( ' ', n, n, u, ldu, uz, lda )
791 *
792  CALL shseqr( 'S', 'V', n, ilo, ihi, t1, lda, wr1, wi1, uz,
793  $ ldu, work, nwork, iinfo )
794  IF( iinfo.NE.0 .AND. iinfo.LE.n+2 ) THEN
795  WRITE( nounit, fmt = 9999 )'SHSEQR(V)', iinfo, n, jtype,
796  $ ioldsd
797  info = abs( iinfo )
798  go to 250
799  END IF
800 *
801 * Compute Z = U' UZ
802 *
803  CALL sgemm( 'T', 'N', n, n, n, one, u, ldu, uz, ldu, zero,
804  $ z, ldu )
805  ntest = 8
806 *
807 * Do Tests 3: | H - Z T Z' | / ( |H| n ulp )
808 * and 4: | I - Z Z' | / ( n ulp )
809 *
810  CALL shst01( n, ilo, ihi, h, lda, t1, lda, z, ldu, work,
811  $ nwork, result( 3 ) )
812 *
813 * Do Tests 5: | A - UZ T (UZ)' | / ( |A| n ulp )
814 * and 6: | I - UZ (UZ)' | / ( n ulp )
815 *
816  CALL shst01( n, ilo, ihi, a, lda, t1, lda, uz, ldu, work,
817  $ nwork, result( 5 ) )
818 *
819 * Do Test 7: | T2 - T1 | / ( |T| n ulp )
820 *
821  CALL sget10( n, n, t2, lda, t1, lda, work, result( 7 ) )
822 *
823 * Do Test 8: | W3 - W1 | / ( max(|W1|,|W3|) ulp )
824 *
825  temp1 = zero
826  temp2 = zero
827  DO 130 j = 1, n
828  temp1 = max( temp1, abs( wr1( j ) )+abs( wi1( j ) ),
829  $ abs( wr3( j ) )+abs( wi3( j ) ) )
830  temp2 = max( temp2, abs( wr1( j )-wr3( j ) )+
831  $ abs( wr1( j )-wr3( j ) ) )
832  130 continue
833 *
834  result( 8 ) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
835 *
836 * Compute the Left and Right Eigenvectors of T
837 *
838 * Compute the Right eigenvector Matrix:
839 *
840  ntest = 9
841  result( 9 ) = ulpinv
842 *
843 * Select last max(N/4,1) real, max(N/4,1) complex eigenvectors
844 *
845  nselc = 0
846  nselr = 0
847  j = n
848  140 continue
849  IF( wi1( j ).EQ.zero ) THEN
850  IF( nselr.LT.max( n / 4, 1 ) ) THEN
851  nselr = nselr + 1
852  SELECT( j ) = .true.
853  ELSE
854  SELECT( j ) = .false.
855  END IF
856  j = j - 1
857  ELSE
858  IF( nselc.LT.max( n / 4, 1 ) ) THEN
859  nselc = nselc + 1
860  SELECT( j ) = .true.
861  SELECT( j-1 ) = .false.
862  ELSE
863  SELECT( j ) = .false.
864  SELECT( j-1 ) = .false.
865  END IF
866  j = j - 2
867  END IF
868  IF( j.GT.0 )
869  $ go to 140
870 *
871  CALL strevc( 'Right', 'All', SELECT, n, t1, lda, dumma, ldu,
872  $ evectr, ldu, n, in, work, iinfo )
873  IF( iinfo.NE.0 ) THEN
874  WRITE( nounit, fmt = 9999 )'STREVC(R,A)', iinfo, n,
875  $ jtype, ioldsd
876  info = abs( iinfo )
877  go to 250
878  END IF
879 *
880 * Test 9: | TR - RW | / ( |T| |R| ulp )
881 *
882  CALL sget22( 'N', 'N', 'N', n, t1, lda, evectr, ldu, wr1,
883  $ wi1, work, dumma( 1 ) )
884  result( 9 ) = dumma( 1 )
885  IF( dumma( 2 ).GT.thresh ) THEN
886  WRITE( nounit, fmt = 9998 )'Right', 'STREVC',
887  $ dumma( 2 ), n, jtype, ioldsd
888  END IF
889 *
890 * Compute selected right eigenvectors and confirm that
891 * they agree with previous right eigenvectors
892 *
893  CALL strevc( 'Right', 'Some', SELECT, n, t1, lda, dumma,
894  $ ldu, evectl, ldu, n, in, work, iinfo )
895  IF( iinfo.NE.0 ) THEN
896  WRITE( nounit, fmt = 9999 )'STREVC(R,S)', iinfo, n,
897  $ jtype, ioldsd
898  info = abs( iinfo )
899  go to 250
900  END IF
901 *
902  k = 1
903  match = .true.
904  DO 170 j = 1, n
905  IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
906  DO 150 jj = 1, n
907  IF( evectr( jj, j ).NE.evectl( jj, k ) ) THEN
908  match = .false.
909  go to 180
910  END IF
911  150 continue
912  k = k + 1
913  ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
914  DO 160 jj = 1, n
915  IF( evectr( jj, j ).NE.evectl( jj, k ) .OR.
916  $ evectr( jj, j+1 ).NE.evectl( jj, k+1 ) ) THEN
917  match = .false.
918  go to 180
919  END IF
920  160 continue
921  k = k + 2
922  END IF
923  170 continue
924  180 continue
925  IF( .NOT.match )
926  $ WRITE( nounit, fmt = 9997 )'Right', 'STREVC', n, jtype,
927  $ ioldsd
928 *
929 * Compute the Left eigenvector Matrix:
930 *
931  ntest = 10
932  result( 10 ) = ulpinv
933  CALL strevc( 'Left', 'All', SELECT, n, t1, lda, evectl, ldu,
934  $ dumma, ldu, n, in, work, iinfo )
935  IF( iinfo.NE.0 ) THEN
936  WRITE( nounit, fmt = 9999 )'STREVC(L,A)', iinfo, n,
937  $ jtype, ioldsd
938  info = abs( iinfo )
939  go to 250
940  END IF
941 *
942 * Test 10: | LT - WL | / ( |T| |L| ulp )
943 *
944  CALL sget22( 'Trans', 'N', 'Conj', n, t1, lda, evectl, ldu,
945  $ wr1, wi1, work, dumma( 3 ) )
946  result( 10 ) = dumma( 3 )
947  IF( dumma( 4 ).GT.thresh ) THEN
948  WRITE( nounit, fmt = 9998 )'Left', 'STREVC', dumma( 4 ),
949  $ n, jtype, ioldsd
950  END IF
951 *
952 * Compute selected left eigenvectors and confirm that
953 * they agree with previous left eigenvectors
954 *
955  CALL strevc( 'Left', 'Some', SELECT, n, t1, lda, evectr,
956  $ ldu, dumma, ldu, n, in, work, iinfo )
957  IF( iinfo.NE.0 ) THEN
958  WRITE( nounit, fmt = 9999 )'STREVC(L,S)', iinfo, n,
959  $ jtype, ioldsd
960  info = abs( iinfo )
961  go to 250
962  END IF
963 *
964  k = 1
965  match = .true.
966  DO 210 j = 1, n
967  IF( SELECT( j ) .AND. wi1( j ).EQ.zero ) THEN
968  DO 190 jj = 1, n
969  IF( evectl( jj, j ).NE.evectr( jj, k ) ) THEN
970  match = .false.
971  go to 220
972  END IF
973  190 continue
974  k = k + 1
975  ELSE IF( SELECT( j ) .AND. wi1( j ).NE.zero ) THEN
976  DO 200 jj = 1, n
977  IF( evectl( jj, j ).NE.evectr( jj, k ) .OR.
978  $ evectl( jj, j+1 ).NE.evectr( jj, k+1 ) ) THEN
979  match = .false.
980  go to 220
981  END IF
982  200 continue
983  k = k + 2
984  END IF
985  210 continue
986  220 continue
987  IF( .NOT.match )
988  $ WRITE( nounit, fmt = 9997 )'Left', 'STREVC', n, jtype,
989  $ ioldsd
990 *
991 * Call SHSEIN for Right eigenvectors of H, do test 11
992 *
993  ntest = 11
994  result( 11 ) = ulpinv
995  DO 230 j = 1, n
996  SELECT( j ) = .true.
997  230 continue
998 *
999  CALL shsein( 'Right', 'Qr', 'Ninitv', SELECT, n, h, lda,
1000  $ wr3, wi3, dumma, ldu, evectx, ldu, n1, in,
1001  $ work, iwork, iwork, iinfo )
1002  IF( iinfo.NE.0 ) THEN
1003  WRITE( nounit, fmt = 9999 )'SHSEIN(R)', iinfo, n, jtype,
1004  $ ioldsd
1005  info = abs( iinfo )
1006  IF( iinfo.LT.0 )
1007  $ go to 250
1008  ELSE
1009 *
1010 * Test 11: | HX - XW | / ( |H| |X| ulp )
1011 *
1012 * (from inverse iteration)
1013 *
1014  CALL sget22( 'N', 'N', 'N', n, h, lda, evectx, ldu, wr3,
1015  $ wi3, work, dumma( 1 ) )
1016  IF( dumma( 1 ).LT.ulpinv )
1017  $ result( 11 ) = dumma( 1 )*aninv
1018  IF( dumma( 2 ).GT.thresh ) THEN
1019  WRITE( nounit, fmt = 9998 )'Right', 'SHSEIN',
1020  $ dumma( 2 ), n, jtype, ioldsd
1021  END IF
1022  END IF
1023 *
1024 * Call SHSEIN for Left eigenvectors of H, do test 12
1025 *
1026  ntest = 12
1027  result( 12 ) = ulpinv
1028  DO 240 j = 1, n
1029  SELECT( j ) = .true.
1030  240 continue
1031 *
1032  CALL shsein( 'Left', 'Qr', 'Ninitv', SELECT, n, h, lda, wr3,
1033  $ wi3, evecty, ldu, dumma, ldu, n1, in, work,
1034  $ iwork, iwork, iinfo )
1035  IF( iinfo.NE.0 ) THEN
1036  WRITE( nounit, fmt = 9999 )'SHSEIN(L)', iinfo, n, jtype,
1037  $ ioldsd
1038  info = abs( iinfo )
1039  IF( iinfo.LT.0 )
1040  $ go to 250
1041  ELSE
1042 *
1043 * Test 12: | YH - WY | / ( |H| |Y| ulp )
1044 *
1045 * (from inverse iteration)
1046 *
1047  CALL sget22( 'C', 'N', 'C', n, h, lda, evecty, ldu, wr3,
1048  $ wi3, work, dumma( 3 ) )
1049  IF( dumma( 3 ).LT.ulpinv )
1050  $ result( 12 ) = dumma( 3 )*aninv
1051  IF( dumma( 4 ).GT.thresh ) THEN
1052  WRITE( nounit, fmt = 9998 )'Left', 'SHSEIN',
1053  $ dumma( 4 ), n, jtype, ioldsd
1054  END IF
1055  END IF
1056 *
1057 * Call SORMHR for Right eigenvectors of A, do test 13
1058 *
1059  ntest = 13
1060  result( 13 ) = ulpinv
1061 *
1062  CALL sormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1063  $ ldu, tau, evectx, ldu, work, nwork, iinfo )
1064  IF( iinfo.NE.0 ) THEN
1065  WRITE( nounit, fmt = 9999 )'SORMHR(R)', iinfo, n, jtype,
1066  $ ioldsd
1067  info = abs( iinfo )
1068  IF( iinfo.LT.0 )
1069  $ go to 250
1070  ELSE
1071 *
1072 * Test 13: | AX - XW | / ( |A| |X| ulp )
1073 *
1074 * (from inverse iteration)
1075 *
1076  CALL sget22( 'N', 'N', 'N', n, a, lda, evectx, ldu, wr3,
1077  $ wi3, work, dumma( 1 ) )
1078  IF( dumma( 1 ).LT.ulpinv )
1079  $ result( 13 ) = dumma( 1 )*aninv
1080  END IF
1081 *
1082 * Call SORMHR for Left eigenvectors of A, do test 14
1083 *
1084  ntest = 14
1085  result( 14 ) = ulpinv
1086 *
1087  CALL sormhr( 'Left', 'No transpose', n, n, ilo, ihi, uu,
1088  $ ldu, tau, evecty, ldu, work, nwork, iinfo )
1089  IF( iinfo.NE.0 ) THEN
1090  WRITE( nounit, fmt = 9999 )'SORMHR(L)', iinfo, n, jtype,
1091  $ ioldsd
1092  info = abs( iinfo )
1093  IF( iinfo.LT.0 )
1094  $ go to 250
1095  ELSE
1096 *
1097 * Test 14: | YA - WY | / ( |A| |Y| ulp )
1098 *
1099 * (from inverse iteration)
1100 *
1101  CALL sget22( 'C', 'N', 'C', n, a, lda, evecty, ldu, wr3,
1102  $ wi3, work, dumma( 3 ) )
1103  IF( dumma( 3 ).LT.ulpinv )
1104  $ result( 14 ) = dumma( 3 )*aninv
1105  END IF
1106 *
1107 * End of Loop -- Check for RESULT(j) > THRESH
1108 *
1109  250 continue
1110 *
1111  ntestt = ntestt + ntest
1112  CALL slafts( 'SHS', n, n, jtype, ntest, result, ioldsd,
1113  $ thresh, nounit, nerrs )
1114 *
1115  260 continue
1116  270 continue
1117 *
1118 * Summary
1119 *
1120  CALL slasum( 'SHS', nounit, nerrs, ntestt )
1121 *
1122  return
1123 *
1124  9999 format( ' SCHKHS: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
1125  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1126  9998 format( ' SCHKHS: ', a, ' Eigenvectors from ', a, ' incorrectly ',
1127  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
1128  $ 'N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5,
1129  $ ')' )
1130  9997 format( ' SCHKHS: Selected ', a, ' Eigenvectors from ', a,
1131  $ ' do not match other eigenvectors ', 9x, 'N=', i6,
1132  $ ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
1133 *
1134 * End of SCHKHS
1135 *
1136  END