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.

zdrgvx.f 25 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757
  1. *> \brief \b ZDRGVX
  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 ZDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
  12. * ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
  13. * S, DTRU, DIF, DIFTRU, WORK, LWORK, RWORK,
  14. * IWORK, LIWORK, RESULT, BWORK, INFO )
  15. *
  16. * .. Scalar Arguments ..
  17. * INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
  18. * $ NSIZE
  19. * DOUBLE PRECISION THRESH
  20. * ..
  21. * .. Array Arguments ..
  22. * LOGICAL BWORK( * )
  23. * INTEGER IWORK( * )
  24. * DOUBLE PRECISION DIF( * ), DIFTRU( * ), DTRU( * ), LSCALE( * ),
  25. * $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * )
  26. * COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ),
  27. * $ B( LDA, * ), BETA( * ), BI( LDA, * ),
  28. * $ VL( LDA, * ), VR( LDA, * ), WORK( * )
  29. * ..
  30. *
  31. *
  32. *> \par Purpose:
  33. * =============
  34. *>
  35. *> \verbatim
  36. *>
  37. *> ZDRGVX checks the nonsymmetric generalized eigenvalue problem
  38. *> expert driver ZGGEVX.
  39. *>
  40. *> ZGGEVX computes the generalized eigenvalues, (optionally) the left
  41. *> and/or right eigenvectors, (optionally) computes a balancing
  42. *> transformation to improve the conditioning, and (optionally)
  43. *> reciprocal condition numbers for the eigenvalues and eigenvectors.
  44. *>
  45. *> When ZDRGVX is called with NSIZE > 0, two types of test matrix pairs
  46. *> are generated by the subroutine DLATM6 and test the driver ZGGEVX.
  47. *> The test matrices have the known exact condition numbers for
  48. *> eigenvalues. For the condition numbers of the eigenvectors
  49. *> corresponding the first and last eigenvalues are also know
  50. *> ``exactly'' (see ZLATM6).
  51. *> For each matrix pair, the following tests will be performed and
  52. *> compared with the threshold THRESH.
  53. *>
  54. *> (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
  55. *>
  56. *> | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
  57. *>
  58. *> where l**H is the conjugate transpose of l.
  59. *>
  60. *> (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
  61. *>
  62. *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
  63. *>
  64. *> (3) The condition number S(i) of eigenvalues computed by ZGGEVX
  65. *> differs less than a factor THRESH from the exact S(i) (see
  66. *> ZLATM6).
  67. *>
  68. *> (4) DIF(i) computed by ZTGSNA differs less than a factor 10*THRESH
  69. *> from the exact value (for the 1st and 5th vectors only).
  70. *>
  71. *> Test Matrices
  72. *> =============
  73. *>
  74. *> Two kinds of test matrix pairs
  75. *> (A, B) = inverse(YH) * (Da, Db) * inverse(X)
  76. *> are used in the tests:
  77. *>
  78. *> 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
  79. *> 0 2+a 0 0 0 0 1 0 0 0
  80. *> 0 0 3+a 0 0 0 0 1 0 0
  81. *> 0 0 0 4+a 0 0 0 0 1 0
  82. *> 0 0 0 0 5+a , 0 0 0 0 1 , and
  83. *>
  84. *> 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0
  85. *> 1 1 0 0 0 0 1 0 0 0
  86. *> 0 0 1 0 0 0 0 1 0 0
  87. *> 0 0 0 1+a 1+b 0 0 0 1 0
  88. *> 0 0 0 -1-b 1+a , 0 0 0 0 1 .
  89. *>
  90. *> In both cases the same inverse(YH) and inverse(X) are used to compute
  91. *> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
  92. *>
  93. *> YH: = 1 0 -y y -y X = 1 0 -x -x x
  94. *> 0 1 -y y -y 0 1 x -x -x
  95. *> 0 0 1 0 0 0 0 1 0 0
  96. *> 0 0 0 1 0 0 0 0 1 0
  97. *> 0 0 0 0 1, 0 0 0 0 1 , where
  98. *>
  99. *> a, b, x and y will have all values independently of each other from
  100. *> { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }.
  101. *> \endverbatim
  102. *
  103. * Arguments:
  104. * ==========
  105. *
  106. *> \param[in] NSIZE
  107. *> \verbatim
  108. *> NSIZE is INTEGER
  109. *> The number of sizes of matrices to use. NSIZE must be at
  110. *> least zero. If it is zero, no randomly generated matrices
  111. *> are tested, but any test matrices read from NIN will be
  112. *> tested. If it is not zero, then N = 5.
  113. *> \endverbatim
  114. *>
  115. *> \param[in] THRESH
  116. *> \verbatim
  117. *> THRESH is DOUBLE PRECISION
  118. *> A test will count as "failed" if the "error", computed as
  119. *> described above, exceeds THRESH. Note that the error
  120. *> is scaled to be O(1), so THRESH should be a reasonably
  121. *> small multiple of 1, e.g., 10 or 100. In particular,
  122. *> it should not depend on the precision (single vs. double)
  123. *> or the size of the matrix. It must be at least zero.
  124. *> \endverbatim
  125. *>
  126. *> \param[in] NIN
  127. *> \verbatim
  128. *> NIN is INTEGER
  129. *> The FORTRAN unit number for reading in the data file of
  130. *> problems to solve.
  131. *> \endverbatim
  132. *>
  133. *> \param[in] NOUT
  134. *> \verbatim
  135. *> NOUT is INTEGER
  136. *> The FORTRAN unit number for printing out error messages
  137. *> (e.g., if a routine returns IINFO not equal to 0.)
  138. *> \endverbatim
  139. *>
  140. *> \param[out] A
  141. *> \verbatim
  142. *> A is COMPLEX*16 array, dimension (LDA, NSIZE)
  143. *> Used to hold the matrix whose eigenvalues are to be
  144. *> computed. On exit, A contains the last matrix actually used.
  145. *> \endverbatim
  146. *>
  147. *> \param[in] LDA
  148. *> \verbatim
  149. *> LDA is INTEGER
  150. *> The leading dimension of A, B, AI, BI, Ao, and Bo.
  151. *> It must be at least 1 and at least NSIZE.
  152. *> \endverbatim
  153. *>
  154. *> \param[out] B
  155. *> \verbatim
  156. *> B is COMPLEX*16 array, dimension (LDA, NSIZE)
  157. *> Used to hold the matrix whose eigenvalues are to be
  158. *> computed. On exit, B contains the last matrix actually used.
  159. *> \endverbatim
  160. *>
  161. *> \param[out] AI
  162. *> \verbatim
  163. *> AI is COMPLEX*16 array, dimension (LDA, NSIZE)
  164. *> Copy of A, modified by ZGGEVX.
  165. *> \endverbatim
  166. *>
  167. *> \param[out] BI
  168. *> \verbatim
  169. *> BI is COMPLEX*16 array, dimension (LDA, NSIZE)
  170. *> Copy of B, modified by ZGGEVX.
  171. *> \endverbatim
  172. *>
  173. *> \param[out] ALPHA
  174. *> \verbatim
  175. *> ALPHA is COMPLEX*16 array, dimension (NSIZE)
  176. *> \endverbatim
  177. *>
  178. *> \param[out] BETA
  179. *> \verbatim
  180. *> BETA is COMPLEX*16 array, dimension (NSIZE)
  181. *>
  182. *> On exit, ALPHA/BETA are the eigenvalues.
  183. *> \endverbatim
  184. *>
  185. *> \param[out] VL
  186. *> \verbatim
  187. *> VL is COMPLEX*16 array, dimension (LDA, NSIZE)
  188. *> VL holds the left eigenvectors computed by ZGGEVX.
  189. *> \endverbatim
  190. *>
  191. *> \param[out] VR
  192. *> \verbatim
  193. *> VR is COMPLEX*16 array, dimension (LDA, NSIZE)
  194. *> VR holds the right eigenvectors computed by ZGGEVX.
  195. *> \endverbatim
  196. *>
  197. *> \param[out] ILO
  198. *> \verbatim
  199. *> ILO is INTEGER
  200. *> \endverbatim
  201. *>
  202. *> \param[out] IHI
  203. *> \verbatim
  204. *> IHI is INTEGER
  205. *> \endverbatim
  206. *>
  207. *> \param[out] LSCALE
  208. *> \verbatim
  209. *> LSCALE is DOUBLE PRECISION array, dimension (N)
  210. *> \endverbatim
  211. *>
  212. *> \param[out] RSCALE
  213. *> \verbatim
  214. *> RSCALE is DOUBLE PRECISION array, dimension (N)
  215. *> \endverbatim
  216. *>
  217. *> \param[out] S
  218. *> \verbatim
  219. *> S is DOUBLE PRECISION array, dimension (N)
  220. *> \endverbatim
  221. *>
  222. *> \param[out] DTRU
  223. *> \verbatim
  224. *> DTRU is DOUBLE PRECISION array, dimension (N)
  225. *> \endverbatim
  226. *>
  227. *> \param[out] DIF
  228. *> \verbatim
  229. *> DIF is DOUBLE PRECISION array, dimension (N)
  230. *> \endverbatim
  231. *>
  232. *> \param[out] DIFTRU
  233. *> \verbatim
  234. *> DIFTRU is DOUBLE PRECISION array, dimension (N)
  235. *> \endverbatim
  236. *>
  237. *> \param[out] WORK
  238. *> \verbatim
  239. *> WORK is COMPLEX*16 array, dimension (LWORK)
  240. *> \endverbatim
  241. *>
  242. *> \param[in] LWORK
  243. *> \verbatim
  244. *> LWORK is INTEGER
  245. *> Leading dimension of WORK. LWORK >= 2*N*N + 2*N
  246. *> \endverbatim
  247. *>
  248. *> \param[out] RWORK
  249. *> \verbatim
  250. *> RWORK is DOUBLE PRECISION array, dimension (6*N)
  251. *> \endverbatim
  252. *>
  253. *> \param[out] IWORK
  254. *> \verbatim
  255. *> IWORK is INTEGER array, dimension (LIWORK)
  256. *> \endverbatim
  257. *>
  258. *> \param[in] LIWORK
  259. *> \verbatim
  260. *> LIWORK is INTEGER
  261. *> Leading dimension of IWORK. LIWORK >= N+2.
  262. *> \endverbatim
  263. *>
  264. *> \param[out] RESULT
  265. *> \verbatim
  266. *> RESULT is DOUBLE PRECISION array, dimension (4)
  267. *> \endverbatim
  268. *>
  269. *> \param[out] BWORK
  270. *> \verbatim
  271. *> BWORK is LOGICAL array, dimension (N)
  272. *> \endverbatim
  273. *>
  274. *> \param[out] INFO
  275. *> \verbatim
  276. *> INFO is INTEGER
  277. *> = 0: successful exit
  278. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  279. *> > 0: A routine returned an error code.
  280. *> \endverbatim
  281. *
  282. * Authors:
  283. * ========
  284. *
  285. *> \author Univ. of Tennessee
  286. *> \author Univ. of California Berkeley
  287. *> \author Univ. of Colorado Denver
  288. *> \author NAG Ltd.
  289. *
  290. *> \ingroup complex16_eig
  291. *
  292. * =====================================================================
  293. SUBROUTINE ZDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
  294. $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
  295. $ S, DTRU, DIF, DIFTRU, WORK, LWORK, RWORK,
  296. $ IWORK, LIWORK, RESULT, BWORK, INFO )
  297. *
  298. * -- LAPACK test routine --
  299. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  300. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  301. *
  302. * .. Scalar Arguments ..
  303. INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
  304. $ NSIZE
  305. DOUBLE PRECISION THRESH
  306. * ..
  307. * .. Array Arguments ..
  308. LOGICAL BWORK( * )
  309. INTEGER IWORK( * )
  310. DOUBLE PRECISION DIF( * ), DIFTRU( * ), DTRU( * ), LSCALE( * ),
  311. $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * )
  312. COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ),
  313. $ B( LDA, * ), BETA( * ), BI( LDA, * ),
  314. $ VL( LDA, * ), VR( LDA, * ), WORK( * )
  315. * ..
  316. *
  317. * =====================================================================
  318. *
  319. * .. Parameters ..
  320. DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
  321. PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1,
  322. $ TNTH = 1.0D-1, HALF = 0.5D+0 )
  323. * ..
  324. * .. Local Scalars ..
  325. INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
  326. $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
  327. DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
  328. $ ULP, ULPINV
  329. * ..
  330. * .. Local Arrays ..
  331. COMPLEX*16 WEIGHT( 5 )
  332. * ..
  333. * .. External Functions ..
  334. INTEGER ILAENV
  335. DOUBLE PRECISION DLAMCH, ZLANGE
  336. EXTERNAL ILAENV, DLAMCH, ZLANGE
  337. * ..
  338. * .. External Subroutines ..
  339. EXTERNAL ALASVM, XERBLA, ZGET52, ZGGEVX, ZLACPY, ZLATM6
  340. * ..
  341. * .. Intrinsic Functions ..
  342. INTRINSIC ABS, DCMPLX, MAX, SQRT
  343. * ..
  344. * .. Executable Statements ..
  345. *
  346. * Check for errors
  347. *
  348. INFO = 0
  349. *
  350. NMAX = 5
  351. *
  352. IF( NSIZE.LT.0 ) THEN
  353. INFO = -1
  354. ELSE IF( THRESH.LT.ZERO ) THEN
  355. INFO = -2
  356. ELSE IF( NIN.LE.0 ) THEN
  357. INFO = -3
  358. ELSE IF( NOUT.LE.0 ) THEN
  359. INFO = -4
  360. ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
  361. INFO = -6
  362. ELSE IF( LIWORK.LT.NMAX+2 ) THEN
  363. INFO = -26
  364. END IF
  365. *
  366. * Compute workspace
  367. * (Note: Comments in the code beginning "Workspace:" describe the
  368. * minimal amount of workspace needed at that point in the code,
  369. * as well as the preferred amount for good performance.
  370. * NB refers to the optimal block size for the immediately
  371. * following subroutine, as returned by ILAENV.)
  372. *
  373. MINWRK = 1
  374. IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
  375. MINWRK = 2*NMAX*( NMAX+1 )
  376. MAXWRK = NMAX*( 1+ILAENV( 1, 'ZGEQRF', ' ', NMAX, 1, NMAX,
  377. $ 0 ) )
  378. MAXWRK = MAX( MAXWRK, 2*NMAX*( NMAX+1 ) )
  379. WORK( 1 ) = MAXWRK
  380. END IF
  381. *
  382. IF( LWORK.LT.MINWRK )
  383. $ INFO = -23
  384. *
  385. IF( INFO.NE.0 ) THEN
  386. CALL XERBLA( 'ZDRGVX', -INFO )
  387. RETURN
  388. END IF
  389. *
  390. N = 5
  391. ULP = DLAMCH( 'P' )
  392. ULPINV = ONE / ULP
  393. THRSH2 = TEN*THRESH
  394. NERRS = 0
  395. NPTKNT = 0
  396. NTESTT = 0
  397. *
  398. IF( NSIZE.EQ.0 )
  399. $ GO TO 90
  400. *
  401. * Parameters used for generating test matrices.
  402. *
  403. WEIGHT( 1 ) = DCMPLX( TNTH, ZERO )
  404. WEIGHT( 2 ) = DCMPLX( HALF, ZERO )
  405. WEIGHT( 3 ) = ONE
  406. WEIGHT( 4 ) = ONE / WEIGHT( 2 )
  407. WEIGHT( 5 ) = ONE / WEIGHT( 1 )
  408. *
  409. DO 80 IPTYPE = 1, 2
  410. DO 70 IWA = 1, 5
  411. DO 60 IWB = 1, 5
  412. DO 50 IWX = 1, 5
  413. DO 40 IWY = 1, 5
  414. *
  415. * generated a pair of test matrix
  416. *
  417. CALL ZLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL,
  418. $ LDA, WEIGHT( IWA ), WEIGHT( IWB ),
  419. $ WEIGHT( IWX ), WEIGHT( IWY ), DTRU,
  420. $ DIFTRU )
  421. *
  422. * Compute eigenvalues/eigenvectors of (A, B).
  423. * Compute eigenvalue/eigenvector condition numbers
  424. * using computed eigenvectors.
  425. *
  426. CALL ZLACPY( 'F', N, N, A, LDA, AI, LDA )
  427. CALL ZLACPY( 'F', N, N, B, LDA, BI, LDA )
  428. *
  429. CALL ZGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI,
  430. $ LDA, ALPHA, BETA, VL, LDA, VR, LDA,
  431. $ ILO, IHI, LSCALE, RSCALE, ANORM,
  432. $ BNORM, S, DIF, WORK, LWORK, RWORK,
  433. $ IWORK, BWORK, LINFO )
  434. IF( LINFO.NE.0 ) THEN
  435. WRITE( NOUT, FMT = 9999 )'ZGGEVX', LINFO, N,
  436. $ IPTYPE, IWA, IWB, IWX, IWY
  437. GO TO 30
  438. END IF
  439. *
  440. * Compute the norm(A, B)
  441. *
  442. CALL ZLACPY( 'Full', N, N, AI, LDA, WORK, N )
  443. CALL ZLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ),
  444. $ N )
  445. ABNORM = ZLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
  446. *
  447. * Tests (1) and (2)
  448. *
  449. RESULT( 1 ) = ZERO
  450. CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA,
  451. $ ALPHA, BETA, WORK, RWORK,
  452. $ RESULT( 1 ) )
  453. IF( RESULT( 2 ).GT.THRESH ) THEN
  454. WRITE( NOUT, FMT = 9998 )'Left', 'ZGGEVX',
  455. $ RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY
  456. END IF
  457. *
  458. RESULT( 2 ) = ZERO
  459. CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA,
  460. $ ALPHA, BETA, WORK, RWORK,
  461. $ RESULT( 2 ) )
  462. IF( RESULT( 3 ).GT.THRESH ) THEN
  463. WRITE( NOUT, FMT = 9998 )'Right', 'ZGGEVX',
  464. $ RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY
  465. END IF
  466. *
  467. * Test (3)
  468. *
  469. RESULT( 3 ) = ZERO
  470. DO 10 I = 1, N
  471. IF( S( I ).EQ.ZERO ) THEN
  472. IF( DTRU( I ).GT.ABNORM*ULP )
  473. $ RESULT( 3 ) = ULPINV
  474. ELSE IF( DTRU( I ).EQ.ZERO ) THEN
  475. IF( S( I ).GT.ABNORM*ULP )
  476. $ RESULT( 3 ) = ULPINV
  477. ELSE
  478. RWORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ),
  479. $ ABS( S( I ) / DTRU( I ) ) )
  480. RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
  481. END IF
  482. 10 CONTINUE
  483. *
  484. * Test (4)
  485. *
  486. RESULT( 4 ) = ZERO
  487. IF( DIF( 1 ).EQ.ZERO ) THEN
  488. IF( DIFTRU( 1 ).GT.ABNORM*ULP )
  489. $ RESULT( 4 ) = ULPINV
  490. ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
  491. IF( DIF( 1 ).GT.ABNORM*ULP )
  492. $ RESULT( 4 ) = ULPINV
  493. ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
  494. IF( DIFTRU( 5 ).GT.ABNORM*ULP )
  495. $ RESULT( 4 ) = ULPINV
  496. ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
  497. IF( DIF( 5 ).GT.ABNORM*ULP )
  498. $ RESULT( 4 ) = ULPINV
  499. ELSE
  500. RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
  501. $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
  502. RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
  503. $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
  504. RESULT( 4 ) = MAX( RATIO1, RATIO2 )
  505. END IF
  506. *
  507. NTESTT = NTESTT + 4
  508. *
  509. * Print out tests which fail.
  510. *
  511. DO 20 J = 1, 4
  512. IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR.
  513. $ ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) )
  514. $ THEN
  515. *
  516. * If this is the first test to fail,
  517. * print a header to the data file.
  518. *
  519. IF( NERRS.EQ.0 ) THEN
  520. WRITE( NOUT, FMT = 9997 )'ZXV'
  521. *
  522. * Print out messages for built-in examples
  523. *
  524. * Matrix types
  525. *
  526. WRITE( NOUT, FMT = 9995 )
  527. WRITE( NOUT, FMT = 9994 )
  528. WRITE( NOUT, FMT = 9993 )
  529. *
  530. * Tests performed
  531. *
  532. WRITE( NOUT, FMT = 9992 )'''',
  533. $ 'transpose', ''''
  534. *
  535. END IF
  536. NERRS = NERRS + 1
  537. IF( RESULT( J ).LT.10000.0D0 ) THEN
  538. WRITE( NOUT, FMT = 9991 )IPTYPE, IWA,
  539. $ IWB, IWX, IWY, J, RESULT( J )
  540. ELSE
  541. WRITE( NOUT, FMT = 9990 )IPTYPE, IWA,
  542. $ IWB, IWX, IWY, J, RESULT( J )
  543. END IF
  544. END IF
  545. 20 CONTINUE
  546. *
  547. 30 CONTINUE
  548. *
  549. 40 CONTINUE
  550. 50 CONTINUE
  551. 60 CONTINUE
  552. 70 CONTINUE
  553. 80 CONTINUE
  554. *
  555. GO TO 150
  556. *
  557. 90 CONTINUE
  558. *
  559. * Read in data from file to check accuracy of condition estimation
  560. * Read input data until N=0
  561. *
  562. READ( NIN, FMT = *, END = 150 )N
  563. IF( N.EQ.0 )
  564. $ GO TO 150
  565. DO 100 I = 1, N
  566. READ( NIN, FMT = * )( A( I, J ), J = 1, N )
  567. 100 CONTINUE
  568. DO 110 I = 1, N
  569. READ( NIN, FMT = * )( B( I, J ), J = 1, N )
  570. 110 CONTINUE
  571. READ( NIN, FMT = * )( DTRU( I ), I = 1, N )
  572. READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N )
  573. *
  574. NPTKNT = NPTKNT + 1
  575. *
  576. * Compute eigenvalues/eigenvectors of (A, B).
  577. * Compute eigenvalue/eigenvector condition numbers
  578. * using computed eigenvectors.
  579. *
  580. CALL ZLACPY( 'F', N, N, A, LDA, AI, LDA )
  581. CALL ZLACPY( 'F', N, N, B, LDA, BI, LDA )
  582. *
  583. CALL ZGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHA, BETA,
  584. $ VL, LDA, VR, LDA, ILO, IHI, LSCALE, RSCALE, ANORM,
  585. $ BNORM, S, DIF, WORK, LWORK, RWORK, IWORK, BWORK,
  586. $ LINFO )
  587. *
  588. IF( LINFO.NE.0 ) THEN
  589. WRITE( NOUT, FMT = 9987 )'ZGGEVX', LINFO, N, NPTKNT
  590. GO TO 140
  591. END IF
  592. *
  593. * Compute the norm(A, B)
  594. *
  595. CALL ZLACPY( 'Full', N, N, AI, LDA, WORK, N )
  596. CALL ZLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N )
  597. ABNORM = ZLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
  598. *
  599. * Tests (1) and (2)
  600. *
  601. RESULT( 1 ) = ZERO
  602. CALL ZGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHA, BETA,
  603. $ WORK, RWORK, RESULT( 1 ) )
  604. IF( RESULT( 2 ).GT.THRESH ) THEN
  605. WRITE( NOUT, FMT = 9986 )'Left', 'ZGGEVX', RESULT( 2 ), N,
  606. $ NPTKNT
  607. END IF
  608. *
  609. RESULT( 2 ) = ZERO
  610. CALL ZGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHA, BETA,
  611. $ WORK, RWORK, RESULT( 2 ) )
  612. IF( RESULT( 3 ).GT.THRESH ) THEN
  613. WRITE( NOUT, FMT = 9986 )'Right', 'ZGGEVX', RESULT( 3 ), N,
  614. $ NPTKNT
  615. END IF
  616. *
  617. * Test (3)
  618. *
  619. RESULT( 3 ) = ZERO
  620. DO 120 I = 1, N
  621. IF( S( I ).EQ.ZERO ) THEN
  622. IF( DTRU( I ).GT.ABNORM*ULP )
  623. $ RESULT( 3 ) = ULPINV
  624. ELSE IF( DTRU( I ).EQ.ZERO ) THEN
  625. IF( S( I ).GT.ABNORM*ULP )
  626. $ RESULT( 3 ) = ULPINV
  627. ELSE
  628. RWORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ),
  629. $ ABS( S( I ) / DTRU( I ) ) )
  630. RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
  631. END IF
  632. 120 CONTINUE
  633. *
  634. * Test (4)
  635. *
  636. RESULT( 4 ) = ZERO
  637. IF( DIF( 1 ).EQ.ZERO ) THEN
  638. IF( DIFTRU( 1 ).GT.ABNORM*ULP )
  639. $ RESULT( 4 ) = ULPINV
  640. ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
  641. IF( DIF( 1 ).GT.ABNORM*ULP )
  642. $ RESULT( 4 ) = ULPINV
  643. ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
  644. IF( DIFTRU( 5 ).GT.ABNORM*ULP )
  645. $ RESULT( 4 ) = ULPINV
  646. ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
  647. IF( DIF( 5 ).GT.ABNORM*ULP )
  648. $ RESULT( 4 ) = ULPINV
  649. ELSE
  650. RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
  651. $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
  652. RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
  653. $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
  654. RESULT( 4 ) = MAX( RATIO1, RATIO2 )
  655. END IF
  656. *
  657. NTESTT = NTESTT + 4
  658. *
  659. * Print out tests which fail.
  660. *
  661. DO 130 J = 1, 4
  662. IF( RESULT( J ).GE.THRSH2 ) THEN
  663. *
  664. * If this is the first test to fail,
  665. * print a header to the data file.
  666. *
  667. IF( NERRS.EQ.0 ) THEN
  668. WRITE( NOUT, FMT = 9997 )'ZXV'
  669. *
  670. * Print out messages for built-in examples
  671. *
  672. * Matrix types
  673. *
  674. WRITE( NOUT, FMT = 9996 )
  675. *
  676. * Tests performed
  677. *
  678. WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
  679. *
  680. END IF
  681. NERRS = NERRS + 1
  682. IF( RESULT( J ).LT.10000.0D0 ) THEN
  683. WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
  684. ELSE
  685. WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
  686. END IF
  687. END IF
  688. 130 CONTINUE
  689. *
  690. 140 CONTINUE
  691. *
  692. GO TO 90
  693. 150 CONTINUE
  694. *
  695. * Summary
  696. *
  697. CALL ALASVM( 'ZXV', NOUT, NERRS, NTESTT, 0 )
  698. *
  699. WORK( 1 ) = MAXWRK
  700. *
  701. RETURN
  702. *
  703. 9999 FORMAT( ' ZDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  704. $ I6, ', JTYPE=', I6, ')' )
  705. *
  706. 9998 FORMAT( ' ZDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
  707. $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
  708. $ 'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5,
  709. $ ', IWX=', I5, ', IWY=', I5 )
  710. *
  711. 9997 FORMAT( / 1X, A3, ' -- Complex Expert Eigenvalue/vector',
  712. $ ' problem driver' )
  713. *
  714. 9996 FORMAT( 'Input Example' )
  715. *
  716. 9995 FORMAT( ' Matrix types: ', / )
  717. *
  718. 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
  719. $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
  720. $ / ' YH and X are left and right eigenvectors. ', / )
  721. *
  722. 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
  723. $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
  724. $ / ' YH and X are left and right eigenvectors. ', / )
  725. *
  726. 9992 FORMAT( / ' Tests performed: ', / 4X,
  727. $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
  728. $ ' r is a right eigenvector and ', A, ' means ', A, '.',
  729. $ / ' 1 = max | ( b A - a B )', A, ' l | / const.',
  730. $ / ' 2 = max | ( b A - a B ) r | / const.',
  731. $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
  732. $ ' over all eigenvalues', /
  733. $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
  734. $ ' over the 1st and 5th eigenvectors', / )
  735. *
  736. 9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
  737. $ I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 )
  738. *
  739. 9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
  740. $ I2, ', IWY=', I2, ', result ', I2, ' is', 1P, D10.3 )
  741. *
  742. 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
  743. $ ' result ', I2, ' is', 0P, F8.2 )
  744. *
  745. 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
  746. $ ' result ', I2, ' is', 1P, D10.3 )
  747. *
  748. 9987 FORMAT( ' ZDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  749. $ I6, ', Input example #', I2, ')' )
  750. *
  751. 9986 FORMAT( ' ZDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
  752. $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
  753. $ 'N=', I6, ', Input Example #', I2, ')' )
  754. *
  755. * End of ZDRGVX
  756. *
  757. END