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.

dchkhs.f 43 kB

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