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.

dget24.f 32 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992
  1. *> \brief \b DGET24
  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 DGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
  12. * H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
  13. * LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
  14. * RESULT, WORK, LWORK, IWORK, BWORK, INFO )
  15. *
  16. * .. Scalar Arguments ..
  17. * LOGICAL COMP
  18. * INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
  19. * DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
  20. * ..
  21. * .. Array Arguments ..
  22. * LOGICAL BWORK( * )
  23. * INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
  24. * DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
  25. * $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
  26. * $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
  27. * $ WR( * ), WRT( * ), WRTMP( * )
  28. * ..
  29. *
  30. *
  31. *> \par Purpose:
  32. * =============
  33. *>
  34. *> \verbatim
  35. *>
  36. *> DGET24 checks the nonsymmetric eigenvalue (Schur form) problem
  37. *> expert driver DGEESX.
  38. *>
  39. *> If COMP = .FALSE., the first 13 of the following tests will be
  40. *> be performed on the input matrix A, and also tests 14 and 15
  41. *> if LWORK is sufficiently large.
  42. *> If COMP = .TRUE., all 17 test will be performed.
  43. *>
  44. *> (1) 0 if T is in Schur form, 1/ulp otherwise
  45. *> (no sorting of eigenvalues)
  46. *>
  47. *> (2) | A - VS T VS' | / ( n |A| ulp )
  48. *>
  49. *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
  50. *> form (no sorting of eigenvalues).
  51. *>
  52. *> (3) | I - VS VS' | / ( n ulp ) (no sorting of eigenvalues).
  53. *>
  54. *> (4) 0 if WR+sqrt(-1)*WI are eigenvalues of T
  55. *> 1/ulp otherwise
  56. *> (no sorting of eigenvalues)
  57. *>
  58. *> (5) 0 if T(with VS) = T(without VS),
  59. *> 1/ulp otherwise
  60. *> (no sorting of eigenvalues)
  61. *>
  62. *> (6) 0 if eigenvalues(with VS) = eigenvalues(without VS),
  63. *> 1/ulp otherwise
  64. *> (no sorting of eigenvalues)
  65. *>
  66. *> (7) 0 if T is in Schur form, 1/ulp otherwise
  67. *> (with sorting of eigenvalues)
  68. *>
  69. *> (8) | A - VS T VS' | / ( n |A| ulp )
  70. *>
  71. *> Here VS is the matrix of Schur eigenvectors, and T is in Schur
  72. *> form (with sorting of eigenvalues).
  73. *>
  74. *> (9) | I - VS VS' | / ( n ulp ) (with sorting of eigenvalues).
  75. *>
  76. *> (10) 0 if WR+sqrt(-1)*WI are eigenvalues of T
  77. *> 1/ulp otherwise
  78. *> If workspace sufficient, also compare WR, WI with and
  79. *> without reciprocal condition numbers
  80. *> (with sorting of eigenvalues)
  81. *>
  82. *> (11) 0 if T(with VS) = T(without VS),
  83. *> 1/ulp otherwise
  84. *> If workspace sufficient, also compare T with and without
  85. *> reciprocal condition numbers
  86. *> (with sorting of eigenvalues)
  87. *>
  88. *> (12) 0 if eigenvalues(with VS) = eigenvalues(without VS),
  89. *> 1/ulp otherwise
  90. *> If workspace sufficient, also compare VS with and without
  91. *> reciprocal condition numbers
  92. *> (with sorting of eigenvalues)
  93. *>
  94. *> (13) if sorting worked and SDIM is the number of
  95. *> eigenvalues which were SELECTed
  96. *> If workspace sufficient, also compare SDIM with and
  97. *> without reciprocal condition numbers
  98. *>
  99. *> (14) if RCONDE the same no matter if VS and/or RCONDV computed
  100. *>
  101. *> (15) if RCONDV the same no matter if VS and/or RCONDE computed
  102. *>
  103. *> (16) |RCONDE - RCDEIN| / cond(RCONDE)
  104. *>
  105. *> RCONDE is the reciprocal average eigenvalue condition number
  106. *> computed by DGEESX and RCDEIN (the precomputed true value)
  107. *> is supplied as input. cond(RCONDE) is the condition number
  108. *> of RCONDE, and takes errors in computing RCONDE into account,
  109. *> so that the resulting quantity should be O(ULP). cond(RCONDE)
  110. *> is essentially given by norm(A)/RCONDV.
  111. *>
  112. *> (17) |RCONDV - RCDVIN| / cond(RCONDV)
  113. *>
  114. *> RCONDV is the reciprocal right invariant subspace condition
  115. *> number computed by DGEESX and RCDVIN (the precomputed true
  116. *> value) is supplied as input. cond(RCONDV) is the condition
  117. *> number of RCONDV, and takes errors in computing RCONDV into
  118. *> account, so that the resulting quantity should be O(ULP).
  119. *> cond(RCONDV) is essentially given by norm(A)/RCONDE.
  120. *> \endverbatim
  121. *
  122. * Arguments:
  123. * ==========
  124. *
  125. *> \param[in] COMP
  126. *> \verbatim
  127. *> COMP is LOGICAL
  128. *> COMP describes which input tests to perform:
  129. *> = .FALSE. if the computed condition numbers are not to
  130. *> be tested against RCDVIN and RCDEIN
  131. *> = .TRUE. if they are to be compared
  132. *> \endverbatim
  133. *>
  134. *> \param[in] JTYPE
  135. *> \verbatim
  136. *> JTYPE is INTEGER
  137. *> Type of input matrix. Used to label output if error occurs.
  138. *> \endverbatim
  139. *>
  140. *> \param[in] ISEED
  141. *> \verbatim
  142. *> ISEED is INTEGER array, dimension (4)
  143. *> If COMP = .FALSE., the random number generator seed
  144. *> used to produce matrix.
  145. *> If COMP = .TRUE., ISEED(1) = the number of the example.
  146. *> Used to label output if error occurs.
  147. *> \endverbatim
  148. *>
  149. *> \param[in] THRESH
  150. *> \verbatim
  151. *> THRESH is DOUBLE PRECISION
  152. *> A test will count as "failed" if the "error", computed as
  153. *> described above, exceeds THRESH. Note that the error
  154. *> is scaled to be O(1), so THRESH should be a reasonably
  155. *> small multiple of 1, e.g., 10 or 100. In particular,
  156. *> it should not depend on the precision (single vs. double)
  157. *> or the size of the matrix. It must be at least zero.
  158. *> \endverbatim
  159. *>
  160. *> \param[in] NOUNIT
  161. *> \verbatim
  162. *> NOUNIT is INTEGER
  163. *> The FORTRAN unit number for printing out error messages
  164. *> (e.g., if a routine returns INFO not equal to 0.)
  165. *> \endverbatim
  166. *>
  167. *> \param[in] N
  168. *> \verbatim
  169. *> N is INTEGER
  170. *> The dimension of A. N must be at least 0.
  171. *> \endverbatim
  172. *>
  173. *> \param[in,out] A
  174. *> \verbatim
  175. *> A is DOUBLE PRECISION array, dimension (LDA, N)
  176. *> Used to hold the matrix whose eigenvalues are to be
  177. *> computed.
  178. *> \endverbatim
  179. *>
  180. *> \param[in] LDA
  181. *> \verbatim
  182. *> LDA is INTEGER
  183. *> The leading dimension of A, and H. LDA must be at
  184. *> least 1 and at least N.
  185. *> \endverbatim
  186. *>
  187. *> \param[out] H
  188. *> \verbatim
  189. *> H is DOUBLE PRECISION array, dimension (LDA, N)
  190. *> Another copy of the test matrix A, modified by DGEESX.
  191. *> \endverbatim
  192. *>
  193. *> \param[out] HT
  194. *> \verbatim
  195. *> HT is DOUBLE PRECISION array, dimension (LDA, N)
  196. *> Yet another copy of the test matrix A, modified by DGEESX.
  197. *> \endverbatim
  198. *>
  199. *> \param[out] WR
  200. *> \verbatim
  201. *> WR is DOUBLE PRECISION array, dimension (N)
  202. *> \endverbatim
  203. *>
  204. *> \param[out] WI
  205. *> \verbatim
  206. *> WI is DOUBLE PRECISION array, dimension (N)
  207. *>
  208. *> The real and imaginary parts of the eigenvalues of A.
  209. *> On exit, WR + WI*i are the eigenvalues of the matrix in A.
  210. *> \endverbatim
  211. *>
  212. *> \param[out] WRT
  213. *> \verbatim
  214. *> WRT is DOUBLE PRECISION array, dimension (N)
  215. *> \endverbatim
  216. *>
  217. *> \param[out] WIT
  218. *> \verbatim
  219. *> WIT is DOUBLE PRECISION array, dimension (N)
  220. *>
  221. *> Like WR, WI, these arrays contain the eigenvalues of A,
  222. *> but those computed when DGEESX only computes a partial
  223. *> eigendecomposition, i.e. not Schur vectors
  224. *> \endverbatim
  225. *>
  226. *> \param[out] WRTMP
  227. *> \verbatim
  228. *> WRTMP is DOUBLE PRECISION array, dimension (N)
  229. *> \endverbatim
  230. *>
  231. *> \param[out] WITMP
  232. *> \verbatim
  233. *> WITMP is DOUBLE PRECISION array, dimension (N)
  234. *>
  235. *> Like WR, WI, these arrays contain the eigenvalues of A,
  236. *> but sorted by increasing real part.
  237. *> \endverbatim
  238. *>
  239. *> \param[out] VS
  240. *> \verbatim
  241. *> VS is DOUBLE PRECISION array, dimension (LDVS, N)
  242. *> VS holds the computed Schur vectors.
  243. *> \endverbatim
  244. *>
  245. *> \param[in] LDVS
  246. *> \verbatim
  247. *> LDVS is INTEGER
  248. *> Leading dimension of VS. Must be at least max(1, N).
  249. *> \endverbatim
  250. *>
  251. *> \param[out] VS1
  252. *> \verbatim
  253. *> VS1 is DOUBLE PRECISION array, dimension (LDVS, N)
  254. *> VS1 holds another copy of the computed Schur vectors.
  255. *> \endverbatim
  256. *>
  257. *> \param[in] RCDEIN
  258. *> \verbatim
  259. *> RCDEIN is DOUBLE PRECISION
  260. *> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
  261. *> condition number for the average of selected eigenvalues.
  262. *> \endverbatim
  263. *>
  264. *> \param[in] RCDVIN
  265. *> \verbatim
  266. *> RCDVIN is DOUBLE PRECISION
  267. *> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
  268. *> condition number for the selected right invariant subspace.
  269. *> \endverbatim
  270. *>
  271. *> \param[in] NSLCT
  272. *> \verbatim
  273. *> NSLCT is INTEGER
  274. *> When COMP = .TRUE. the number of selected eigenvalues
  275. *> corresponding to the precomputed values RCDEIN and RCDVIN.
  276. *> \endverbatim
  277. *>
  278. *> \param[in] ISLCT
  279. *> \verbatim
  280. *> ISLCT is INTEGER array, dimension (NSLCT)
  281. *> When COMP = .TRUE. ISLCT selects the eigenvalues of the
  282. *> input matrix corresponding to the precomputed values RCDEIN
  283. *> and RCDVIN. For I=1, ... ,NSLCT, if ISLCT(I) = J, then the
  284. *> eigenvalue with the J-th largest real part is selected.
  285. *> Not referenced if COMP = .FALSE.
  286. *> \endverbatim
  287. *>
  288. *> \param[out] RESULT
  289. *> \verbatim
  290. *> RESULT is DOUBLE PRECISION array, dimension (17)
  291. *> The values computed by the 17 tests described above.
  292. *> The values are currently limited to 1/ulp, to avoid
  293. *> overflow.
  294. *> \endverbatim
  295. *>
  296. *> \param[out] WORK
  297. *> \verbatim
  298. *> WORK is DOUBLE PRECISION array, dimension (LWORK)
  299. *> \endverbatim
  300. *>
  301. *> \param[in] LWORK
  302. *> \verbatim
  303. *> LWORK is INTEGER
  304. *> The number of entries in WORK to be passed to DGEESX. This
  305. *> must be at least 3*N, and N+N**2 if tests 14--16 are to
  306. *> be performed.
  307. *> \endverbatim
  308. *>
  309. *> \param[out] IWORK
  310. *> \verbatim
  311. *> IWORK is INTEGER array, dimension (N*N)
  312. *> \endverbatim
  313. *>
  314. *> \param[out] BWORK
  315. *> \verbatim
  316. *> BWORK is LOGICAL array, dimension (N)
  317. *> \endverbatim
  318. *>
  319. *> \param[out] INFO
  320. *> \verbatim
  321. *> INFO is INTEGER
  322. *> If 0, successful exit.
  323. *> If <0, input parameter -INFO had an incorrect value.
  324. *> If >0, DGEESX returned an error code, the absolute
  325. *> value of which is returned.
  326. *> \endverbatim
  327. *
  328. * Authors:
  329. * ========
  330. *
  331. *> \author Univ. of Tennessee
  332. *> \author Univ. of California Berkeley
  333. *> \author Univ. of Colorado Denver
  334. *> \author NAG Ltd.
  335. *
  336. *> \ingroup double_eig
  337. *
  338. * =====================================================================
  339. SUBROUTINE DGET24( COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA,
  340. $ H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS,
  341. $ LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT,
  342. $ RESULT, WORK, LWORK, IWORK, BWORK, INFO )
  343. *
  344. * -- LAPACK test routine --
  345. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  346. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  347. *
  348. * .. Scalar Arguments ..
  349. LOGICAL COMP
  350. INTEGER INFO, JTYPE, LDA, LDVS, LWORK, N, NOUNIT, NSLCT
  351. DOUBLE PRECISION RCDEIN, RCDVIN, THRESH
  352. * ..
  353. * .. Array Arguments ..
  354. LOGICAL BWORK( * )
  355. INTEGER ISEED( 4 ), ISLCT( * ), IWORK( * )
  356. DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
  357. $ RESULT( 17 ), VS( LDVS, * ), VS1( LDVS, * ),
  358. $ WI( * ), WIT( * ), WITMP( * ), WORK( * ),
  359. $ WR( * ), WRT( * ), WRTMP( * )
  360. * ..
  361. *
  362. * =====================================================================
  363. *
  364. * .. Parameters ..
  365. DOUBLE PRECISION ZERO, ONE
  366. PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
  367. DOUBLE PRECISION EPSIN
  368. PARAMETER ( EPSIN = 5.9605D-8 )
  369. * ..
  370. * .. Local Scalars ..
  371. CHARACTER SORT
  372. INTEGER I, IINFO, ISORT, ITMP, J, KMIN, KNTEIG, LIWORK,
  373. $ RSUB, SDIM, SDIM1
  374. DOUBLE PRECISION ANORM, EPS, RCNDE1, RCNDV1, RCONDE, RCONDV,
  375. $ SMLNUM, TMP, TOL, TOLIN, ULP, ULPINV, V, VIMIN,
  376. $ VRMIN, WNORM
  377. * ..
  378. * .. Local Arrays ..
  379. INTEGER IPNT( 20 )
  380. * ..
  381. * .. Arrays in Common ..
  382. LOGICAL SELVAL( 20 )
  383. DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
  384. * ..
  385. * .. Scalars in Common ..
  386. INTEGER SELDIM, SELOPT
  387. * ..
  388. * .. Common blocks ..
  389. COMMON / SSLCT / SELOPT, SELDIM, SELVAL, SELWR, SELWI
  390. * ..
  391. * .. External Functions ..
  392. LOGICAL DSLECT
  393. DOUBLE PRECISION DLAMCH, DLANGE
  394. EXTERNAL DSLECT, DLAMCH, DLANGE
  395. * ..
  396. * .. External Subroutines ..
  397. EXTERNAL DCOPY, DGEESX, DGEMM, DLACPY, DORT01, XERBLA
  398. * ..
  399. * .. Intrinsic Functions ..
  400. INTRINSIC ABS, DBLE, MAX, MIN, SIGN, SQRT
  401. * ..
  402. * .. Executable Statements ..
  403. *
  404. * Check for errors
  405. *
  406. INFO = 0
  407. IF( THRESH.LT.ZERO ) THEN
  408. INFO = -3
  409. ELSE IF( NOUNIT.LE.0 ) THEN
  410. INFO = -5
  411. ELSE IF( N.LT.0 ) THEN
  412. INFO = -6
  413. ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
  414. INFO = -8
  415. ELSE IF( LDVS.LT.1 .OR. LDVS.LT.N ) THEN
  416. INFO = -18
  417. ELSE IF( LWORK.LT.3*N ) THEN
  418. INFO = -26
  419. END IF
  420. *
  421. IF( INFO.NE.0 ) THEN
  422. CALL XERBLA( 'DGET24', -INFO )
  423. RETURN
  424. END IF
  425. *
  426. * Quick return if nothing to do
  427. *
  428. DO 10 I = 1, 17
  429. RESULT( I ) = -ONE
  430. 10 CONTINUE
  431. *
  432. IF( N.EQ.0 )
  433. $ RETURN
  434. *
  435. * Important constants
  436. *
  437. SMLNUM = DLAMCH( 'Safe minimum' )
  438. ULP = DLAMCH( 'Precision' )
  439. ULPINV = ONE / ULP
  440. *
  441. * Perform tests (1)-(13)
  442. *
  443. SELOPT = 0
  444. LIWORK = N*N
  445. DO 120 ISORT = 0, 1
  446. IF( ISORT.EQ.0 ) THEN
  447. SORT = 'N'
  448. RSUB = 0
  449. ELSE
  450. SORT = 'S'
  451. RSUB = 6
  452. END IF
  453. *
  454. * Compute Schur form and Schur vectors, and test them
  455. *
  456. CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
  457. CALL DGEESX( 'V', SORT, DSLECT, 'N', N, H, LDA, SDIM, WR, WI,
  458. $ VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
  459. $ LIWORK, BWORK, IINFO )
  460. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  461. RESULT( 1+RSUB ) = ULPINV
  462. IF( JTYPE.NE.22 ) THEN
  463. WRITE( NOUNIT, FMT = 9998 )'DGEESX1', IINFO, N, JTYPE,
  464. $ ISEED
  465. ELSE
  466. WRITE( NOUNIT, FMT = 9999 )'DGEESX1', IINFO, N,
  467. $ ISEED( 1 )
  468. END IF
  469. INFO = ABS( IINFO )
  470. RETURN
  471. END IF
  472. IF( ISORT.EQ.0 ) THEN
  473. CALL DCOPY( N, WR, 1, WRTMP, 1 )
  474. CALL DCOPY( N, WI, 1, WITMP, 1 )
  475. END IF
  476. *
  477. * Do Test (1) or Test (7)
  478. *
  479. RESULT( 1+RSUB ) = ZERO
  480. DO 30 J = 1, N - 2
  481. DO 20 I = J + 2, N
  482. IF( H( I, J ).NE.ZERO )
  483. $ RESULT( 1+RSUB ) = ULPINV
  484. 20 CONTINUE
  485. 30 CONTINUE
  486. DO 40 I = 1, N - 2
  487. IF( H( I+1, I ).NE.ZERO .AND. H( I+2, I+1 ).NE.ZERO )
  488. $ RESULT( 1+RSUB ) = ULPINV
  489. 40 CONTINUE
  490. DO 50 I = 1, N - 1
  491. IF( H( I+1, I ).NE.ZERO ) THEN
  492. IF( H( I, I ).NE.H( I+1, I+1 ) .OR. H( I, I+1 ).EQ.
  493. $ ZERO .OR. SIGN( ONE, H( I+1, I ) ).EQ.
  494. $ SIGN( ONE, H( I, I+1 ) ) )RESULT( 1+RSUB ) = ULPINV
  495. END IF
  496. 50 CONTINUE
  497. *
  498. * Test (2) or (8): Compute norm(A - Q*H*Q') / (norm(A) * N * ULP)
  499. *
  500. * Copy A to VS1, used as workspace
  501. *
  502. CALL DLACPY( ' ', N, N, A, LDA, VS1, LDVS )
  503. *
  504. * Compute Q*H and store in HT.
  505. *
  506. CALL DGEMM( 'No transpose', 'No transpose', N, N, N, ONE, VS,
  507. $ LDVS, H, LDA, ZERO, HT, LDA )
  508. *
  509. * Compute A - Q*H*Q'
  510. *
  511. CALL DGEMM( 'No transpose', 'Transpose', N, N, N, -ONE, HT,
  512. $ LDA, VS, LDVS, ONE, VS1, LDVS )
  513. *
  514. ANORM = MAX( DLANGE( '1', N, N, A, LDA, WORK ), SMLNUM )
  515. WNORM = DLANGE( '1', N, N, VS1, LDVS, WORK )
  516. *
  517. IF( ANORM.GT.WNORM ) THEN
  518. RESULT( 2+RSUB ) = ( WNORM / ANORM ) / ( N*ULP )
  519. ELSE
  520. IF( ANORM.LT.ONE ) THEN
  521. RESULT( 2+RSUB ) = ( MIN( WNORM, N*ANORM ) / ANORM ) /
  522. $ ( N*ULP )
  523. ELSE
  524. RESULT( 2+RSUB ) = MIN( WNORM / ANORM, DBLE( N ) ) /
  525. $ ( N*ULP )
  526. END IF
  527. END IF
  528. *
  529. * Test (3) or (9): Compute norm( I - Q'*Q ) / ( N * ULP )
  530. *
  531. CALL DORT01( 'Columns', N, N, VS, LDVS, WORK, LWORK,
  532. $ RESULT( 3+RSUB ) )
  533. *
  534. * Do Test (4) or Test (10)
  535. *
  536. RESULT( 4+RSUB ) = ZERO
  537. DO 60 I = 1, N
  538. IF( H( I, I ).NE.WR( I ) )
  539. $ RESULT( 4+RSUB ) = ULPINV
  540. 60 CONTINUE
  541. IF( N.GT.1 ) THEN
  542. IF( H( 2, 1 ).EQ.ZERO .AND. WI( 1 ).NE.ZERO )
  543. $ RESULT( 4+RSUB ) = ULPINV
  544. IF( H( N, N-1 ).EQ.ZERO .AND. WI( N ).NE.ZERO )
  545. $ RESULT( 4+RSUB ) = ULPINV
  546. END IF
  547. DO 70 I = 1, N - 1
  548. IF( H( I+1, I ).NE.ZERO ) THEN
  549. TMP = SQRT( ABS( H( I+1, I ) ) )*
  550. $ SQRT( ABS( H( I, I+1 ) ) )
  551. RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
  552. $ ABS( WI( I )-TMP ) /
  553. $ MAX( ULP*TMP, SMLNUM ) )
  554. RESULT( 4+RSUB ) = MAX( RESULT( 4+RSUB ),
  555. $ ABS( WI( I+1 )+TMP ) /
  556. $ MAX( ULP*TMP, SMLNUM ) )
  557. ELSE IF( I.GT.1 ) THEN
  558. IF( H( I+1, I ).EQ.ZERO .AND. H( I, I-1 ).EQ.ZERO .AND.
  559. $ WI( I ).NE.ZERO )RESULT( 4+RSUB ) = ULPINV
  560. END IF
  561. 70 CONTINUE
  562. *
  563. * Do Test (5) or Test (11)
  564. *
  565. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  566. CALL DGEESX( 'N', SORT, DSLECT, 'N', N, HT, LDA, SDIM, WRT,
  567. $ WIT, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK,
  568. $ LIWORK, BWORK, IINFO )
  569. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  570. RESULT( 5+RSUB ) = ULPINV
  571. IF( JTYPE.NE.22 ) THEN
  572. WRITE( NOUNIT, FMT = 9998 )'DGEESX2', IINFO, N, JTYPE,
  573. $ ISEED
  574. ELSE
  575. WRITE( NOUNIT, FMT = 9999 )'DGEESX2', IINFO, N,
  576. $ ISEED( 1 )
  577. END IF
  578. INFO = ABS( IINFO )
  579. GO TO 250
  580. END IF
  581. *
  582. RESULT( 5+RSUB ) = ZERO
  583. DO 90 J = 1, N
  584. DO 80 I = 1, N
  585. IF( H( I, J ).NE.HT( I, J ) )
  586. $ RESULT( 5+RSUB ) = ULPINV
  587. 80 CONTINUE
  588. 90 CONTINUE
  589. *
  590. * Do Test (6) or Test (12)
  591. *
  592. RESULT( 6+RSUB ) = ZERO
  593. DO 100 I = 1, N
  594. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  595. $ RESULT( 6+RSUB ) = ULPINV
  596. 100 CONTINUE
  597. *
  598. * Do Test (13)
  599. *
  600. IF( ISORT.EQ.1 ) THEN
  601. RESULT( 13 ) = ZERO
  602. KNTEIG = 0
  603. DO 110 I = 1, N
  604. IF( DSLECT( WR( I ), WI( I ) ) .OR.
  605. $ DSLECT( WR( I ), -WI( I ) ) )KNTEIG = KNTEIG + 1
  606. IF( I.LT.N ) THEN
  607. IF( ( DSLECT( WR( I+1 ), WI( I+1 ) ) .OR.
  608. $ DSLECT( WR( I+1 ), -WI( I+1 ) ) ) .AND.
  609. $ ( .NOT.( DSLECT( WR( I ),
  610. $ WI( I ) ) .OR. DSLECT( WR( I ),
  611. $ -WI( I ) ) ) ) .AND. IINFO.NE.N+2 )RESULT( 13 )
  612. $ = ULPINV
  613. END IF
  614. 110 CONTINUE
  615. IF( SDIM.NE.KNTEIG )
  616. $ RESULT( 13 ) = ULPINV
  617. END IF
  618. *
  619. 120 CONTINUE
  620. *
  621. * If there is enough workspace, perform tests (14) and (15)
  622. * as well as (10) through (13)
  623. *
  624. IF( LWORK.GE.N+( N*N ) / 2 ) THEN
  625. *
  626. * Compute both RCONDE and RCONDV with VS
  627. *
  628. SORT = 'S'
  629. RESULT( 14 ) = ZERO
  630. RESULT( 15 ) = ZERO
  631. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  632. CALL DGEESX( 'V', SORT, DSLECT, 'B', N, HT, LDA, SDIM1, WRT,
  633. $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
  634. $ IWORK, LIWORK, BWORK, IINFO )
  635. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  636. RESULT( 14 ) = ULPINV
  637. RESULT( 15 ) = ULPINV
  638. IF( JTYPE.NE.22 ) THEN
  639. WRITE( NOUNIT, FMT = 9998 )'DGEESX3', IINFO, N, JTYPE,
  640. $ ISEED
  641. ELSE
  642. WRITE( NOUNIT, FMT = 9999 )'DGEESX3', IINFO, N,
  643. $ ISEED( 1 )
  644. END IF
  645. INFO = ABS( IINFO )
  646. GO TO 250
  647. END IF
  648. *
  649. * Perform tests (10), (11), (12), and (13)
  650. *
  651. DO 140 I = 1, N
  652. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  653. $ RESULT( 10 ) = ULPINV
  654. DO 130 J = 1, N
  655. IF( H( I, J ).NE.HT( I, J ) )
  656. $ RESULT( 11 ) = ULPINV
  657. IF( VS( I, J ).NE.VS1( I, J ) )
  658. $ RESULT( 12 ) = ULPINV
  659. 130 CONTINUE
  660. 140 CONTINUE
  661. IF( SDIM.NE.SDIM1 )
  662. $ RESULT( 13 ) = ULPINV
  663. *
  664. * Compute both RCONDE and RCONDV without VS, and compare
  665. *
  666. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  667. CALL DGEESX( 'N', SORT, DSLECT, 'B', N, HT, LDA, SDIM1, WRT,
  668. $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
  669. $ IWORK, LIWORK, BWORK, IINFO )
  670. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  671. RESULT( 14 ) = ULPINV
  672. RESULT( 15 ) = ULPINV
  673. IF( JTYPE.NE.22 ) THEN
  674. WRITE( NOUNIT, FMT = 9998 )'DGEESX4', IINFO, N, JTYPE,
  675. $ ISEED
  676. ELSE
  677. WRITE( NOUNIT, FMT = 9999 )'DGEESX4', IINFO, N,
  678. $ ISEED( 1 )
  679. END IF
  680. INFO = ABS( IINFO )
  681. GO TO 250
  682. END IF
  683. *
  684. * Perform tests (14) and (15)
  685. *
  686. IF( RCNDE1.NE.RCONDE )
  687. $ RESULT( 14 ) = ULPINV
  688. IF( RCNDV1.NE.RCONDV )
  689. $ RESULT( 15 ) = ULPINV
  690. *
  691. * Perform tests (10), (11), (12), and (13)
  692. *
  693. DO 160 I = 1, N
  694. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  695. $ RESULT( 10 ) = ULPINV
  696. DO 150 J = 1, N
  697. IF( H( I, J ).NE.HT( I, J ) )
  698. $ RESULT( 11 ) = ULPINV
  699. IF( VS( I, J ).NE.VS1( I, J ) )
  700. $ RESULT( 12 ) = ULPINV
  701. 150 CONTINUE
  702. 160 CONTINUE
  703. IF( SDIM.NE.SDIM1 )
  704. $ RESULT( 13 ) = ULPINV
  705. *
  706. * Compute RCONDE with VS, and compare
  707. *
  708. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  709. CALL DGEESX( 'V', SORT, DSLECT, 'E', N, HT, LDA, SDIM1, WRT,
  710. $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
  711. $ IWORK, LIWORK, BWORK, IINFO )
  712. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  713. RESULT( 14 ) = ULPINV
  714. IF( JTYPE.NE.22 ) THEN
  715. WRITE( NOUNIT, FMT = 9998 )'DGEESX5', IINFO, N, JTYPE,
  716. $ ISEED
  717. ELSE
  718. WRITE( NOUNIT, FMT = 9999 )'DGEESX5', IINFO, N,
  719. $ ISEED( 1 )
  720. END IF
  721. INFO = ABS( IINFO )
  722. GO TO 250
  723. END IF
  724. *
  725. * Perform test (14)
  726. *
  727. IF( RCNDE1.NE.RCONDE )
  728. $ RESULT( 14 ) = ULPINV
  729. *
  730. * Perform tests (10), (11), (12), and (13)
  731. *
  732. DO 180 I = 1, N
  733. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  734. $ RESULT( 10 ) = ULPINV
  735. DO 170 J = 1, N
  736. IF( H( I, J ).NE.HT( I, J ) )
  737. $ RESULT( 11 ) = ULPINV
  738. IF( VS( I, J ).NE.VS1( I, J ) )
  739. $ RESULT( 12 ) = ULPINV
  740. 170 CONTINUE
  741. 180 CONTINUE
  742. IF( SDIM.NE.SDIM1 )
  743. $ RESULT( 13 ) = ULPINV
  744. *
  745. * Compute RCONDE without VS, and compare
  746. *
  747. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  748. CALL DGEESX( 'N', SORT, DSLECT, 'E', N, HT, LDA, SDIM1, WRT,
  749. $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
  750. $ IWORK, LIWORK, BWORK, IINFO )
  751. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  752. RESULT( 14 ) = ULPINV
  753. IF( JTYPE.NE.22 ) THEN
  754. WRITE( NOUNIT, FMT = 9998 )'DGEESX6', IINFO, N, JTYPE,
  755. $ ISEED
  756. ELSE
  757. WRITE( NOUNIT, FMT = 9999 )'DGEESX6', IINFO, N,
  758. $ ISEED( 1 )
  759. END IF
  760. INFO = ABS( IINFO )
  761. GO TO 250
  762. END IF
  763. *
  764. * Perform test (14)
  765. *
  766. IF( RCNDE1.NE.RCONDE )
  767. $ RESULT( 14 ) = ULPINV
  768. *
  769. * Perform tests (10), (11), (12), and (13)
  770. *
  771. DO 200 I = 1, N
  772. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  773. $ RESULT( 10 ) = ULPINV
  774. DO 190 J = 1, N
  775. IF( H( I, J ).NE.HT( I, J ) )
  776. $ RESULT( 11 ) = ULPINV
  777. IF( VS( I, J ).NE.VS1( I, J ) )
  778. $ RESULT( 12 ) = ULPINV
  779. 190 CONTINUE
  780. 200 CONTINUE
  781. IF( SDIM.NE.SDIM1 )
  782. $ RESULT( 13 ) = ULPINV
  783. *
  784. * Compute RCONDV with VS, and compare
  785. *
  786. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  787. CALL DGEESX( 'V', SORT, DSLECT, 'V', N, HT, LDA, SDIM1, WRT,
  788. $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
  789. $ IWORK, LIWORK, BWORK, IINFO )
  790. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  791. RESULT( 15 ) = ULPINV
  792. IF( JTYPE.NE.22 ) THEN
  793. WRITE( NOUNIT, FMT = 9998 )'DGEESX7', IINFO, N, JTYPE,
  794. $ ISEED
  795. ELSE
  796. WRITE( NOUNIT, FMT = 9999 )'DGEESX7', IINFO, N,
  797. $ ISEED( 1 )
  798. END IF
  799. INFO = ABS( IINFO )
  800. GO TO 250
  801. END IF
  802. *
  803. * Perform test (15)
  804. *
  805. IF( RCNDV1.NE.RCONDV )
  806. $ RESULT( 15 ) = ULPINV
  807. *
  808. * Perform tests (10), (11), (12), and (13)
  809. *
  810. DO 220 I = 1, N
  811. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  812. $ RESULT( 10 ) = ULPINV
  813. DO 210 J = 1, N
  814. IF( H( I, J ).NE.HT( I, J ) )
  815. $ RESULT( 11 ) = ULPINV
  816. IF( VS( I, J ).NE.VS1( I, J ) )
  817. $ RESULT( 12 ) = ULPINV
  818. 210 CONTINUE
  819. 220 CONTINUE
  820. IF( SDIM.NE.SDIM1 )
  821. $ RESULT( 13 ) = ULPINV
  822. *
  823. * Compute RCONDV without VS, and compare
  824. *
  825. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  826. CALL DGEESX( 'N', SORT, DSLECT, 'V', N, HT, LDA, SDIM1, WRT,
  827. $ WIT, VS1, LDVS, RCNDE1, RCNDV1, WORK, LWORK,
  828. $ IWORK, LIWORK, BWORK, IINFO )
  829. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  830. RESULT( 15 ) = ULPINV
  831. IF( JTYPE.NE.22 ) THEN
  832. WRITE( NOUNIT, FMT = 9998 )'DGEESX8', IINFO, N, JTYPE,
  833. $ ISEED
  834. ELSE
  835. WRITE( NOUNIT, FMT = 9999 )'DGEESX8', IINFO, N,
  836. $ ISEED( 1 )
  837. END IF
  838. INFO = ABS( IINFO )
  839. GO TO 250
  840. END IF
  841. *
  842. * Perform test (15)
  843. *
  844. IF( RCNDV1.NE.RCONDV )
  845. $ RESULT( 15 ) = ULPINV
  846. *
  847. * Perform tests (10), (11), (12), and (13)
  848. *
  849. DO 240 I = 1, N
  850. IF( WR( I ).NE.WRT( I ) .OR. WI( I ).NE.WIT( I ) )
  851. $ RESULT( 10 ) = ULPINV
  852. DO 230 J = 1, N
  853. IF( H( I, J ).NE.HT( I, J ) )
  854. $ RESULT( 11 ) = ULPINV
  855. IF( VS( I, J ).NE.VS1( I, J ) )
  856. $ RESULT( 12 ) = ULPINV
  857. 230 CONTINUE
  858. 240 CONTINUE
  859. IF( SDIM.NE.SDIM1 )
  860. $ RESULT( 13 ) = ULPINV
  861. *
  862. END IF
  863. *
  864. 250 CONTINUE
  865. *
  866. * If there are precomputed reciprocal condition numbers, compare
  867. * computed values with them.
  868. *
  869. IF( COMP ) THEN
  870. *
  871. * First set up SELOPT, SELDIM, SELVAL, SELWR, and SELWI so that
  872. * the logical function DSLECT selects the eigenvalues specified
  873. * by NSLCT and ISLCT.
  874. *
  875. SELDIM = N
  876. SELOPT = 1
  877. EPS = MAX( ULP, EPSIN )
  878. DO 260 I = 1, N
  879. IPNT( I ) = I
  880. SELVAL( I ) = .FALSE.
  881. SELWR( I ) = WRTMP( I )
  882. SELWI( I ) = WITMP( I )
  883. 260 CONTINUE
  884. DO 280 I = 1, N - 1
  885. KMIN = I
  886. VRMIN = WRTMP( I )
  887. VIMIN = WITMP( I )
  888. DO 270 J = I + 1, N
  889. IF( WRTMP( J ).LT.VRMIN ) THEN
  890. KMIN = J
  891. VRMIN = WRTMP( J )
  892. VIMIN = WITMP( J )
  893. END IF
  894. 270 CONTINUE
  895. WRTMP( KMIN ) = WRTMP( I )
  896. WITMP( KMIN ) = WITMP( I )
  897. WRTMP( I ) = VRMIN
  898. WITMP( I ) = VIMIN
  899. ITMP = IPNT( I )
  900. IPNT( I ) = IPNT( KMIN )
  901. IPNT( KMIN ) = ITMP
  902. 280 CONTINUE
  903. DO 290 I = 1, NSLCT
  904. SELVAL( IPNT( ISLCT( I ) ) ) = .TRUE.
  905. 290 CONTINUE
  906. *
  907. * Compute condition numbers
  908. *
  909. CALL DLACPY( 'F', N, N, A, LDA, HT, LDA )
  910. CALL DGEESX( 'N', 'S', DSLECT, 'B', N, HT, LDA, SDIM1, WRT,
  911. $ WIT, VS1, LDVS, RCONDE, RCONDV, WORK, LWORK,
  912. $ IWORK, LIWORK, BWORK, IINFO )
  913. IF( IINFO.NE.0 .AND. IINFO.NE.N+2 ) THEN
  914. RESULT( 16 ) = ULPINV
  915. RESULT( 17 ) = ULPINV
  916. WRITE( NOUNIT, FMT = 9999 )'DGEESX9', IINFO, N, ISEED( 1 )
  917. INFO = ABS( IINFO )
  918. GO TO 300
  919. END IF
  920. *
  921. * Compare condition number for average of selected eigenvalues
  922. * taking its condition number into account
  923. *
  924. ANORM = DLANGE( '1', N, N, A, LDA, WORK )
  925. V = MAX( DBLE( N )*EPS*ANORM, SMLNUM )
  926. IF( ANORM.EQ.ZERO )
  927. $ V = ONE
  928. IF( V.GT.RCONDV ) THEN
  929. TOL = ONE
  930. ELSE
  931. TOL = V / RCONDV
  932. END IF
  933. IF( V.GT.RCDVIN ) THEN
  934. TOLIN = ONE
  935. ELSE
  936. TOLIN = V / RCDVIN
  937. END IF
  938. TOL = MAX( TOL, SMLNUM / EPS )
  939. TOLIN = MAX( TOLIN, SMLNUM / EPS )
  940. IF( EPS*( RCDEIN-TOLIN ).GT.RCONDE+TOL ) THEN
  941. RESULT( 16 ) = ULPINV
  942. ELSE IF( RCDEIN-TOLIN.GT.RCONDE+TOL ) THEN
  943. RESULT( 16 ) = ( RCDEIN-TOLIN ) / ( RCONDE+TOL )
  944. ELSE IF( RCDEIN+TOLIN.LT.EPS*( RCONDE-TOL ) ) THEN
  945. RESULT( 16 ) = ULPINV
  946. ELSE IF( RCDEIN+TOLIN.LT.RCONDE-TOL ) THEN
  947. RESULT( 16 ) = ( RCONDE-TOL ) / ( RCDEIN+TOLIN )
  948. ELSE
  949. RESULT( 16 ) = ONE
  950. END IF
  951. *
  952. * Compare condition numbers for right invariant subspace
  953. * taking its condition number into account
  954. *
  955. IF( V.GT.RCONDV*RCONDE ) THEN
  956. TOL = RCONDV
  957. ELSE
  958. TOL = V / RCONDE
  959. END IF
  960. IF( V.GT.RCDVIN*RCDEIN ) THEN
  961. TOLIN = RCDVIN
  962. ELSE
  963. TOLIN = V / RCDEIN
  964. END IF
  965. TOL = MAX( TOL, SMLNUM / EPS )
  966. TOLIN = MAX( TOLIN, SMLNUM / EPS )
  967. IF( EPS*( RCDVIN-TOLIN ).GT.RCONDV+TOL ) THEN
  968. RESULT( 17 ) = ULPINV
  969. ELSE IF( RCDVIN-TOLIN.GT.RCONDV+TOL ) THEN
  970. RESULT( 17 ) = ( RCDVIN-TOLIN ) / ( RCONDV+TOL )
  971. ELSE IF( RCDVIN+TOLIN.LT.EPS*( RCONDV-TOL ) ) THEN
  972. RESULT( 17 ) = ULPINV
  973. ELSE IF( RCDVIN+TOLIN.LT.RCONDV-TOL ) THEN
  974. RESULT( 17 ) = ( RCONDV-TOL ) / ( RCDVIN+TOLIN )
  975. ELSE
  976. RESULT( 17 ) = ONE
  977. END IF
  978. *
  979. 300 CONTINUE
  980. *
  981. END IF
  982. *
  983. 9999 FORMAT( ' DGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  984. $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
  985. 9998 FORMAT( ' DGET24: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  986. $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
  987. *
  988. RETURN
  989. *
  990. * End of DGET24
  991. *
  992. END