You can not select more than 25 topics Topics must start with a chinese character,a letter or number, can include dashes ('-') and can be up to 35 characters long.

schkhs.f 43 kB

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