PageRenderTime 52ms CodeModel.GetById 19ms RepoModel.GetById 1ms app.codeStats 0ms

/lapack-netlib/TESTING/EIG/ddrgvx.f

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