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.

sdrvvx.f 36 kB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019
  1. *> \brief \b SDRVVX
  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 SDRVVX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  12. * NIUNIT, NOUNIT, A, LDA, H, WR, WI, WR1, WI1,
  13. * VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV, RCNDV1,
  14. * RCDVIN, RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1,
  15. * RESULT, WORK, NWORK, IWORK, INFO )
  16. *
  17. * .. Scalar Arguments ..
  18. * INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NIUNIT, NOUNIT,
  19. * $ NSIZES, NTYPES, NWORK
  20. * REAL THRESH
  21. * ..
  22. * .. Array Arguments ..
  23. * LOGICAL DOTYPE( * )
  24. * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
  25. * REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
  26. * $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
  27. * $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
  28. * $ RESULT( 11 ), SCALE( * ), SCALE1( * ),
  29. * $ VL( LDVL, * ), VR( LDVR, * ), WI( * ),
  30. * $ WI1( * ), WORK( * ), WR( * ), WR1( * )
  31. * ..
  32. *
  33. *
  34. *> \par Purpose:
  35. * =============
  36. *>
  37. *> \verbatim
  38. *>
  39. *> SDRVVX checks the nonsymmetric eigenvalue problem expert driver
  40. *> SGEEVX.
  41. *>
  42. *> SDRVVX uses both test matrices generated randomly depending on
  43. *> data supplied in the calling sequence, as well as on data
  44. *> read from an input file and including precomputed condition
  45. *> numbers to which it compares the ones it computes.
  46. *>
  47. *> When SDRVVX is called, a number of matrix "sizes" ("n's") and a
  48. *> number of matrix "types" are specified in the calling sequence.
  49. *> For each size ("n") and each type of matrix, one matrix will be
  50. *> generated and used to test the nonsymmetric eigenroutines. For
  51. *> each matrix, 9 tests will be performed:
  52. *>
  53. *> (1) | A * VR - VR * W | / ( n |A| ulp )
  54. *>
  55. *> Here VR is the matrix of unit right eigenvectors.
  56. *> W is a block diagonal matrix, with a 1x1 block for each
  57. *> real eigenvalue and a 2x2 block for each complex conjugate
  58. *> pair. If eigenvalues j and j+1 are a complex conjugate pair,
  59. *> so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
  60. *> 2 x 2 block corresponding to the pair will be:
  61. *>
  62. *> ( wr wi )
  63. *> ( -wi wr )
  64. *>
  65. *> Such a block multiplying an n x 2 matrix ( ur ui ) on the
  66. *> right will be the same as multiplying ur + i*ui by wr + i*wi.
  67. *>
  68. *> (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
  69. *>
  70. *> Here VL is the matrix of unit left eigenvectors, A**H is the
  71. *> conjugate transpose of A, and W is as above.
  72. *>
  73. *> (3) | |VR(i)| - 1 | / ulp and largest component real
  74. *>
  75. *> VR(i) denotes the i-th column of VR.
  76. *>
  77. *> (4) | |VL(i)| - 1 | / ulp and largest component real
  78. *>
  79. *> VL(i) denotes the i-th column of VL.
  80. *>
  81. *> (5) W(full) = W(partial)
  82. *>
  83. *> W(full) denotes the eigenvalues computed when VR, VL, RCONDV
  84. *> and RCONDE are also computed, and W(partial) denotes the
  85. *> eigenvalues computed when only some of VR, VL, RCONDV, and
  86. *> RCONDE are computed.
  87. *>
  88. *> (6) VR(full) = VR(partial)
  89. *>
  90. *> VR(full) denotes the right eigenvectors computed when VL, RCONDV
  91. *> and RCONDE are computed, and VR(partial) denotes the result
  92. *> when only some of VL and RCONDV are computed.
  93. *>
  94. *> (7) VL(full) = VL(partial)
  95. *>
  96. *> VL(full) denotes the left eigenvectors computed when VR, RCONDV
  97. *> and RCONDE are computed, and VL(partial) denotes the result
  98. *> when only some of VR and RCONDV are computed.
  99. *>
  100. *> (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
  101. *> SCALE, ILO, IHI, ABNRM (partial)
  102. *> 1/ulp otherwise
  103. *>
  104. *> SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
  105. *> (full) is when VR, VL, RCONDE and RCONDV are also computed, and
  106. *> (partial) is when some are not computed.
  107. *>
  108. *> (9) RCONDV(full) = RCONDV(partial)
  109. *>
  110. *> RCONDV(full) denotes the reciprocal condition numbers of the
  111. *> right eigenvectors computed when VR, VL and RCONDE are also
  112. *> computed. RCONDV(partial) denotes the reciprocal condition
  113. *> numbers when only some of VR, VL and RCONDE are computed.
  114. *>
  115. *> The "sizes" are specified by an array NN(1:NSIZES); the value of
  116. *> each element NN(j) specifies one size.
  117. *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
  118. *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
  119. *> Currently, the list of possible types is:
  120. *>
  121. *> (1) The zero matrix.
  122. *> (2) The identity matrix.
  123. *> (3) A (transposed) Jordan block, with 1's on the diagonal.
  124. *>
  125. *> (4) A diagonal matrix with evenly spaced entries
  126. *> 1, ..., ULP and random signs.
  127. *> (ULP = (first number larger than 1) - 1 )
  128. *> (5) A diagonal matrix with geometrically spaced entries
  129. *> 1, ..., ULP and random signs.
  130. *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
  131. *> and random signs.
  132. *>
  133. *> (7) Same as (4), but multiplied by a constant near
  134. *> the overflow threshold
  135. *> (8) Same as (4), but multiplied by a constant near
  136. *> the underflow threshold
  137. *>
  138. *> (9) A matrix of the form U' T U, where U is orthogonal and
  139. *> T has evenly spaced entries 1, ..., ULP with random signs
  140. *> on the diagonal and random O(1) entries in the upper
  141. *> triangle.
  142. *>
  143. *> (10) A matrix of the form U' T U, where U is orthogonal and
  144. *> T has geometrically spaced entries 1, ..., ULP with random
  145. *> signs on the diagonal and random O(1) entries in the upper
  146. *> triangle.
  147. *>
  148. *> (11) A matrix of the form U' T U, where U is orthogonal and
  149. *> T has "clustered" entries 1, ULP,..., ULP with random
  150. *> signs on the diagonal and random O(1) entries in the upper
  151. *> triangle.
  152. *>
  153. *> (12) A matrix of the form U' T U, where U is orthogonal and
  154. *> T has real or complex conjugate paired eigenvalues randomly
  155. *> chosen from ( ULP, 1 ) and random O(1) entries in the upper
  156. *> triangle.
  157. *>
  158. *> (13) A matrix of the form X' T X, where X has condition
  159. *> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
  160. *> with random signs on the diagonal and random O(1) entries
  161. *> in the upper triangle.
  162. *>
  163. *> (14) A matrix of the form X' T X, where X has condition
  164. *> SQRT( ULP ) and T has geometrically spaced entries
  165. *> 1, ..., ULP with random signs on the diagonal and random
  166. *> O(1) entries in the upper triangle.
  167. *>
  168. *> (15) A matrix of the form X' T X, where X has condition
  169. *> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
  170. *> with random signs on the diagonal and random O(1) entries
  171. *> in the upper triangle.
  172. *>
  173. *> (16) A matrix of the form X' T X, where X has condition
  174. *> SQRT( ULP ) and T has real or complex conjugate paired
  175. *> eigenvalues randomly chosen from ( ULP, 1 ) and random
  176. *> O(1) entries in the upper triangle.
  177. *>
  178. *> (17) Same as (16), but multiplied by a constant
  179. *> near the overflow threshold
  180. *> (18) Same as (16), but multiplied by a constant
  181. *> near the underflow threshold
  182. *>
  183. *> (19) Nonsymmetric matrix with random entries chosen from (-1,1).
  184. *> If N is at least 4, all entries in first two rows and last
  185. *> row, and first column and last two columns are zero.
  186. *> (20) Same as (19), but multiplied by a constant
  187. *> near the overflow threshold
  188. *> (21) Same as (19), but multiplied by a constant
  189. *> near the underflow threshold
  190. *>
  191. *> In addition, an input file will be read from logical unit number
  192. *> NIUNIT. The file contains matrices along with precomputed
  193. *> eigenvalues and reciprocal condition numbers for the eigenvalues
  194. *> and right eigenvectors. For these matrices, in addition to tests
  195. *> (1) to (9) we will compute the following two tests:
  196. *>
  197. *> (10) |RCONDV - RCDVIN| / cond(RCONDV)
  198. *>
  199. *> RCONDV is the reciprocal right eigenvector condition number
  200. *> computed by SGEEVX and RCDVIN (the precomputed true value)
  201. *> is supplied as input. cond(RCONDV) is the condition number of
  202. *> RCONDV, and takes errors in computing RCONDV into account, so
  203. *> that the resulting quantity should be O(ULP). cond(RCONDV) is
  204. *> essentially given by norm(A)/RCONDE.
  205. *>
  206. *> (11) |RCONDE - RCDEIN| / cond(RCONDE)
  207. *>
  208. *> RCONDE is the reciprocal eigenvalue condition number
  209. *> computed by SGEEVX and RCDEIN (the precomputed true value)
  210. *> is supplied as input. cond(RCONDE) is the condition number
  211. *> of RCONDE, and takes errors in computing RCONDE into account,
  212. *> so that the resulting quantity should be O(ULP). cond(RCONDE)
  213. *> is essentially given by norm(A)/RCONDV.
  214. *> \endverbatim
  215. *
  216. * Arguments:
  217. * ==========
  218. *
  219. *> \param[in] NSIZES
  220. *> \verbatim
  221. *> NSIZES is INTEGER
  222. *> The number of sizes of matrices to use. NSIZES must be at
  223. *> least zero. If it is zero, no randomly generated matrices
  224. *> are tested, but any test matrices read from NIUNIT will be
  225. *> tested.
  226. *> \endverbatim
  227. *>
  228. *> \param[in] NN
  229. *> \verbatim
  230. *> NN is INTEGER array, dimension (NSIZES)
  231. *> An array containing the sizes to be used for the matrices.
  232. *> Zero values will be skipped. The values must be at least
  233. *> zero.
  234. *> \endverbatim
  235. *>
  236. *> \param[in] NTYPES
  237. *> \verbatim
  238. *> NTYPES is INTEGER
  239. *> The number of elements in DOTYPE. NTYPES must be at least
  240. *> zero. If it is zero, no randomly generated test matrices
  241. *> are tested, but and test matrices read from NIUNIT will be
  242. *> tested. If it is MAXTYP+1 and NSIZES is 1, then an
  243. *> additional type, MAXTYP+1 is defined, which is to use
  244. *> whatever matrix is in A. This is only useful if
  245. *> DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
  246. *> \endverbatim
  247. *>
  248. *> \param[in] DOTYPE
  249. *> \verbatim
  250. *> DOTYPE is LOGICAL array, dimension (NTYPES)
  251. *> If DOTYPE(j) is .TRUE., then for each size in NN a
  252. *> matrix of that size and of type j will be generated.
  253. *> If NTYPES is smaller than the maximum number of types
  254. *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
  255. *> MAXTYP will not be generated. If NTYPES is larger
  256. *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
  257. *> will be ignored.
  258. *> \endverbatim
  259. *>
  260. *> \param[in,out] ISEED
  261. *> \verbatim
  262. *> ISEED is INTEGER array, dimension (4)
  263. *> On entry ISEED specifies the seed of the random number
  264. *> generator. The array elements should be between 0 and 4095;
  265. *> if not they will be reduced mod 4096. Also, ISEED(4) must
  266. *> be odd. The random number generator uses a linear
  267. *> congruential sequence limited to small integers, and so
  268. *> should produce machine independent random numbers. The
  269. *> values of ISEED are changed on exit, and can be used in the
  270. *> next call to SDRVVX to continue the same random number
  271. *> sequence.
  272. *> \endverbatim
  273. *>
  274. *> \param[in] THRESH
  275. *> \verbatim
  276. *> THRESH is REAL
  277. *> A test will count as "failed" if the "error", computed as
  278. *> described above, exceeds THRESH. Note that the error
  279. *> is scaled to be O(1), so THRESH should be a reasonably
  280. *> small multiple of 1, e.g., 10 or 100. In particular,
  281. *> it should not depend on the precision (single vs. double)
  282. *> or the size of the matrix. It must be at least zero.
  283. *> \endverbatim
  284. *>
  285. *> \param[in] NIUNIT
  286. *> \verbatim
  287. *> NIUNIT is INTEGER
  288. *> The FORTRAN unit number for reading in the data file of
  289. *> problems to solve.
  290. *> \endverbatim
  291. *>
  292. *> \param[in] NOUNIT
  293. *> \verbatim
  294. *> NOUNIT is INTEGER
  295. *> The FORTRAN unit number for printing out error messages
  296. *> (e.g., if a routine returns INFO not equal to 0.)
  297. *> \endverbatim
  298. *>
  299. *> \param[out] A
  300. *> \verbatim
  301. *> A is REAL array, dimension
  302. *> (LDA, max(NN,12))
  303. *> Used to hold the matrix whose eigenvalues are to be
  304. *> computed. On exit, A contains the last matrix actually used.
  305. *> \endverbatim
  306. *>
  307. *> \param[in] LDA
  308. *> \verbatim
  309. *> LDA is INTEGER
  310. *> The leading dimension of the arrays A and H.
  311. *> LDA >= max(NN,12), since 12 is the dimension of the largest
  312. *> matrix in the precomputed input file.
  313. *> \endverbatim
  314. *>
  315. *> \param[out] H
  316. *> \verbatim
  317. *> H is REAL array, dimension
  318. *> (LDA, max(NN,12))
  319. *> Another copy of the test matrix A, modified by SGEEVX.
  320. *> \endverbatim
  321. *>
  322. *> \param[out] WR
  323. *> \verbatim
  324. *> WR is REAL array, dimension (max(NN))
  325. *> \endverbatim
  326. *>
  327. *> \param[out] WI
  328. *> \verbatim
  329. *> WI is REAL array, dimension (max(NN))
  330. *> The real and imaginary parts of the eigenvalues of A.
  331. *> On exit, WR + WI*i are the eigenvalues of the matrix in A.
  332. *> \endverbatim
  333. *>
  334. *> \param[out] WR1
  335. *> \verbatim
  336. *> WR1 is REAL array, dimension (max(NN,12))
  337. *> \endverbatim
  338. *>
  339. *> \param[out] WI1
  340. *> \verbatim
  341. *> WI1 is REAL array, dimension (max(NN,12))
  342. *>
  343. *> Like WR, WI, these arrays contain the eigenvalues of A,
  344. *> but those computed when SGEEVX only computes a partial
  345. *> eigendecomposition, i.e. not the eigenvalues and left
  346. *> and right eigenvectors.
  347. *> \endverbatim
  348. *>
  349. *> \param[out] VL
  350. *> \verbatim
  351. *> VL is REAL array, dimension
  352. *> (LDVL, max(NN,12))
  353. *> VL holds the computed left eigenvectors.
  354. *> \endverbatim
  355. *>
  356. *> \param[in] LDVL
  357. *> \verbatim
  358. *> LDVL is INTEGER
  359. *> Leading dimension of VL. Must be at least max(1,max(NN,12)).
  360. *> \endverbatim
  361. *>
  362. *> \param[out] VR
  363. *> \verbatim
  364. *> VR is REAL array, dimension
  365. *> (LDVR, max(NN,12))
  366. *> VR holds the computed right eigenvectors.
  367. *> \endverbatim
  368. *>
  369. *> \param[in] LDVR
  370. *> \verbatim
  371. *> LDVR is INTEGER
  372. *> Leading dimension of VR. Must be at least max(1,max(NN,12)).
  373. *> \endverbatim
  374. *>
  375. *> \param[out] LRE
  376. *> \verbatim
  377. *> LRE is REAL array, dimension
  378. *> (LDLRE, max(NN,12))
  379. *> LRE holds the computed right or left eigenvectors.
  380. *> \endverbatim
  381. *>
  382. *> \param[in] LDLRE
  383. *> \verbatim
  384. *> LDLRE is INTEGER
  385. *> Leading dimension of LRE. Must be at least max(1,max(NN,12))
  386. *> \endverbatim
  387. *>
  388. *> \param[out] RCONDV
  389. *> \verbatim
  390. *> RCONDV is REAL array, dimension (N)
  391. *> RCONDV holds the computed reciprocal condition numbers
  392. *> for eigenvectors.
  393. *> \endverbatim
  394. *>
  395. *> \param[out] RCNDV1
  396. *> \verbatim
  397. *> RCNDV1 is REAL array, dimension (N)
  398. *> RCNDV1 holds more computed reciprocal condition numbers
  399. *> for eigenvectors.
  400. *> \endverbatim
  401. *>
  402. *> \param[out] RCDVIN
  403. *> \verbatim
  404. *> RCDVIN is REAL array, dimension (N)
  405. *> When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
  406. *> condition numbers for eigenvectors to be compared with
  407. *> RCONDV.
  408. *> \endverbatim
  409. *>
  410. *> \param[out] RCONDE
  411. *> \verbatim
  412. *> RCONDE is REAL array, dimension (N)
  413. *> RCONDE holds the computed reciprocal condition numbers
  414. *> for eigenvalues.
  415. *> \endverbatim
  416. *>
  417. *> \param[out] RCNDE1
  418. *> \verbatim
  419. *> RCNDE1 is REAL array, dimension (N)
  420. *> RCNDE1 holds more computed reciprocal condition numbers
  421. *> for eigenvalues.
  422. *> \endverbatim
  423. *>
  424. *> \param[out] RCDEIN
  425. *> \verbatim
  426. *> RCDEIN is REAL array, dimension (N)
  427. *> When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
  428. *> condition numbers for eigenvalues to be compared with
  429. *> RCONDE.
  430. *> \endverbatim
  431. *>
  432. *> \param[out] SCALE
  433. *> \verbatim
  434. *> SCALE is REAL array, dimension (N)
  435. *> Holds information describing balancing of matrix.
  436. *> \endverbatim
  437. *>
  438. *> \param[out] SCALE1
  439. *> \verbatim
  440. *> SCALE1 is REAL array, dimension (N)
  441. *> Holds information describing balancing of matrix.
  442. *> \endverbatim
  443. *>
  444. *> \param[out] RESULT
  445. *> \verbatim
  446. *> RESULT is REAL array, dimension (11)
  447. *> The values computed by the seven tests described above.
  448. *> The values are currently limited to 1/ulp, to avoid overflow.
  449. *> \endverbatim
  450. *>
  451. *> \param[out] WORK
  452. *> \verbatim
  453. *> WORK is REAL array, dimension (NWORK)
  454. *> \endverbatim
  455. *>
  456. *> \param[in] NWORK
  457. *> \verbatim
  458. *> NWORK is INTEGER
  459. *> The number of entries in WORK. This must be at least
  460. *> max(6*12+2*12**2,6*NN(j)+2*NN(j)**2) =
  461. *> max( 360 ,6*NN(j)+2*NN(j)**2) for all j.
  462. *> \endverbatim
  463. *>
  464. *> \param[out] IWORK
  465. *> \verbatim
  466. *> IWORK is INTEGER array, dimension (2*max(NN,12))
  467. *> \endverbatim
  468. *>
  469. *> \param[out] INFO
  470. *> \verbatim
  471. *> INFO is INTEGER
  472. *> If 0, then successful exit.
  473. *> If <0, then input parameter -INFO is incorrect.
  474. *> If >0, SLATMR, SLATMS, SLATME or SGET23 returned an error
  475. *> code, and INFO is its absolute value.
  476. *>
  477. *>-----------------------------------------------------------------------
  478. *>
  479. *> Some Local Variables and Parameters:
  480. *> ---- ----- --------- --- ----------
  481. *>
  482. *> ZERO, ONE Real 0 and 1.
  483. *> MAXTYP The number of types defined.
  484. *> NMAX Largest value in NN or 12.
  485. *> NERRS The number of tests which have exceeded THRESH
  486. *> COND, CONDS,
  487. *> IMODE Values to be passed to the matrix generators.
  488. *> ANORM Norm of A; passed to matrix generators.
  489. *>
  490. *> OVFL, UNFL Overflow and underflow thresholds.
  491. *> ULP, ULPINV Finest relative precision and its inverse.
  492. *> RTULP, RTULPI Square roots of the previous 4 values.
  493. *>
  494. *> The following four arrays decode JTYPE:
  495. *> KTYPE(j) The general type (1-10) for type "j".
  496. *> KMODE(j) The MODE value to be passed to the matrix
  497. *> generator for type "j".
  498. *> KMAGN(j) The order of magnitude ( O(1),
  499. *> O(overflow^(1/2) ), O(underflow^(1/2) )
  500. *> KCONDS(j) Selectw whether CONDS is to be 1 or
  501. *> 1/sqrt(ulp). (0 means irrelevant.)
  502. *> \endverbatim
  503. *
  504. * Authors:
  505. * ========
  506. *
  507. *> \author Univ. of Tennessee
  508. *> \author Univ. of California Berkeley
  509. *> \author Univ. of Colorado Denver
  510. *> \author NAG Ltd.
  511. *
  512. *> \date June 2016
  513. *
  514. *> \ingroup single_eig
  515. *
  516. * =====================================================================
  517. SUBROUTINE SDRVVX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  518. $ NIUNIT, NOUNIT, A, LDA, H, WR, WI, WR1, WI1,
  519. $ VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV, RCNDV1,
  520. $ RCDVIN, RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1,
  521. $ RESULT, WORK, NWORK, IWORK, INFO )
  522. *
  523. * -- LAPACK test routine (version 3.7.0) --
  524. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  525. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  526. * June 2016
  527. *
  528. * .. Scalar Arguments ..
  529. INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NIUNIT, NOUNIT,
  530. $ NSIZES, NTYPES, NWORK
  531. REAL THRESH
  532. * ..
  533. * .. Array Arguments ..
  534. LOGICAL DOTYPE( * )
  535. INTEGER ISEED( 4 ), IWORK( * ), NN( * )
  536. REAL A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
  537. $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
  538. $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
  539. $ RESULT( 11 ), SCALE( * ), SCALE1( * ),
  540. $ VL( LDVL, * ), VR( LDVR, * ), WI( * ),
  541. $ WI1( * ), WORK( * ), WR( * ), WR1( * )
  542. * ..
  543. *
  544. * =====================================================================
  545. *
  546. * .. Parameters ..
  547. REAL ZERO, ONE
  548. PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
  549. INTEGER MAXTYP
  550. PARAMETER ( MAXTYP = 21 )
  551. * ..
  552. * .. Local Scalars ..
  553. LOGICAL BADNN
  554. CHARACTER BALANC
  555. CHARACTER*3 PATH
  556. INTEGER I, IBAL, IINFO, IMODE, ITYPE, IWK, J, JCOL,
  557. $ JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL,
  558. $ NMAX, NNWORK, NTEST, NTESTF, NTESTT
  559. REAL ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP,
  560. $ ULPINV, UNFL
  561. * ..
  562. * .. Local Arrays ..
  563. CHARACTER ADUMMA( 1 ), BAL( 4 )
  564. INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
  565. $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
  566. $ KTYPE( MAXTYP )
  567. * ..
  568. * .. External Functions ..
  569. REAL SLAMCH
  570. EXTERNAL SLAMCH
  571. * ..
  572. * .. External Subroutines ..
  573. EXTERNAL SGET23, SLABAD, SLASUM, SLATME, SLATMR, SLATMS,
  574. $ SLASET, XERBLA
  575. * ..
  576. * .. Intrinsic Functions ..
  577. INTRINSIC ABS, MAX, MIN, SQRT
  578. * ..
  579. * .. Data statements ..
  580. DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
  581. DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
  582. $ 3, 1, 2, 3 /
  583. DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
  584. $ 1, 5, 5, 5, 4, 3, 1 /
  585. DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
  586. DATA BAL / 'N', 'P', 'S', 'B' /
  587. * ..
  588. * .. Executable Statements ..
  589. *
  590. PATH( 1: 1 ) = 'Single precision'
  591. PATH( 2: 3 ) = 'VX'
  592. *
  593. * Check for errors
  594. *
  595. NTESTT = 0
  596. NTESTF = 0
  597. INFO = 0
  598. *
  599. * Important constants
  600. *
  601. BADNN = .FALSE.
  602. *
  603. * 12 is the largest dimension in the input file of precomputed
  604. * problems
  605. *
  606. NMAX = 12
  607. DO 10 J = 1, NSIZES
  608. NMAX = MAX( NMAX, NN( J ) )
  609. IF( NN( J ).LT.0 )
  610. $ BADNN = .TRUE.
  611. 10 CONTINUE
  612. *
  613. * Check for errors
  614. *
  615. IF( NSIZES.LT.0 ) THEN
  616. INFO = -1
  617. ELSE IF( BADNN ) THEN
  618. INFO = -2
  619. ELSE IF( NTYPES.LT.0 ) THEN
  620. INFO = -3
  621. ELSE IF( THRESH.LT.ZERO ) THEN
  622. INFO = -6
  623. ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
  624. INFO = -10
  625. ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN
  626. INFO = -17
  627. ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN
  628. INFO = -19
  629. ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN
  630. INFO = -21
  631. ELSE IF( 6*NMAX+2*NMAX**2.GT.NWORK ) THEN
  632. INFO = -32
  633. END IF
  634. *
  635. IF( INFO.NE.0 ) THEN
  636. CALL XERBLA( 'SDRVVX', -INFO )
  637. RETURN
  638. END IF
  639. *
  640. * If nothing to do check on NIUNIT
  641. *
  642. IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
  643. $ GO TO 160
  644. *
  645. * More Important constants
  646. *
  647. UNFL = SLAMCH( 'Safe minimum' )
  648. OVFL = ONE / UNFL
  649. CALL SLABAD( UNFL, OVFL )
  650. ULP = SLAMCH( 'Precision' )
  651. ULPINV = ONE / ULP
  652. RTULP = SQRT( ULP )
  653. RTULPI = ONE / RTULP
  654. *
  655. * Loop over sizes, types
  656. *
  657. NERRS = 0
  658. *
  659. DO 150 JSIZE = 1, NSIZES
  660. N = NN( JSIZE )
  661. IF( NSIZES.NE.1 ) THEN
  662. MTYPES = MIN( MAXTYP, NTYPES )
  663. ELSE
  664. MTYPES = MIN( MAXTYP+1, NTYPES )
  665. END IF
  666. *
  667. DO 140 JTYPE = 1, MTYPES
  668. IF( .NOT.DOTYPE( JTYPE ) )
  669. $ GO TO 140
  670. *
  671. * Save ISEED in case of an error.
  672. *
  673. DO 20 J = 1, 4
  674. IOLDSD( J ) = ISEED( J )
  675. 20 CONTINUE
  676. *
  677. * Compute "A"
  678. *
  679. * Control parameters:
  680. *
  681. * KMAGN KCONDS KMODE KTYPE
  682. * =1 O(1) 1 clustered 1 zero
  683. * =2 large large clustered 2 identity
  684. * =3 small exponential Jordan
  685. * =4 arithmetic diagonal, (w/ eigenvalues)
  686. * =5 random log symmetric, w/ eigenvalues
  687. * =6 random general, w/ eigenvalues
  688. * =7 random diagonal
  689. * =8 random symmetric
  690. * =9 random general
  691. * =10 random triangular
  692. *
  693. IF( MTYPES.GT.MAXTYP )
  694. $ GO TO 90
  695. *
  696. ITYPE = KTYPE( JTYPE )
  697. IMODE = KMODE( JTYPE )
  698. *
  699. * Compute norm
  700. *
  701. GO TO ( 30, 40, 50 )KMAGN( JTYPE )
  702. *
  703. 30 CONTINUE
  704. ANORM = ONE
  705. GO TO 60
  706. *
  707. 40 CONTINUE
  708. ANORM = OVFL*ULP
  709. GO TO 60
  710. *
  711. 50 CONTINUE
  712. ANORM = UNFL*ULPINV
  713. GO TO 60
  714. *
  715. 60 CONTINUE
  716. *
  717. CALL SLASET( 'Full', LDA, N, ZERO, ZERO, A, LDA )
  718. IINFO = 0
  719. COND = ULPINV
  720. *
  721. * Special Matrices -- Identity & Jordan block
  722. *
  723. * Zero
  724. *
  725. IF( ITYPE.EQ.1 ) THEN
  726. IINFO = 0
  727. *
  728. ELSE IF( ITYPE.EQ.2 ) THEN
  729. *
  730. * Identity
  731. *
  732. DO 70 JCOL = 1, N
  733. A( JCOL, JCOL ) = ANORM
  734. 70 CONTINUE
  735. *
  736. ELSE IF( ITYPE.EQ.3 ) THEN
  737. *
  738. * Jordan Block
  739. *
  740. DO 80 JCOL = 1, N
  741. A( JCOL, JCOL ) = ANORM
  742. IF( JCOL.GT.1 )
  743. $ A( JCOL, JCOL-1 ) = ONE
  744. 80 CONTINUE
  745. *
  746. ELSE IF( ITYPE.EQ.4 ) THEN
  747. *
  748. * Diagonal Matrix, [Eigen]values Specified
  749. *
  750. CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
  751. $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
  752. $ IINFO )
  753. *
  754. ELSE IF( ITYPE.EQ.5 ) THEN
  755. *
  756. * Symmetric, eigenvalues specified
  757. *
  758. CALL SLATMS( N, N, 'S', ISEED, 'S', WORK, IMODE, COND,
  759. $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
  760. $ IINFO )
  761. *
  762. ELSE IF( ITYPE.EQ.6 ) THEN
  763. *
  764. * General, eigenvalues specified
  765. *
  766. IF( KCONDS( JTYPE ).EQ.1 ) THEN
  767. CONDS = ONE
  768. ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
  769. CONDS = RTULPI
  770. ELSE
  771. CONDS = ZERO
  772. END IF
  773. *
  774. ADUMMA( 1 ) = ' '
  775. CALL SLATME( N, 'S', ISEED, WORK, IMODE, COND, ONE,
  776. $ ADUMMA, 'T', 'T', 'T', WORK( N+1 ), 4,
  777. $ CONDS, N, N, ANORM, A, LDA, WORK( 2*N+1 ),
  778. $ IINFO )
  779. *
  780. ELSE IF( ITYPE.EQ.7 ) THEN
  781. *
  782. * Diagonal, random eigenvalues
  783. *
  784. CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
  785. $ 'T', 'N', WORK( N+1 ), 1, ONE,
  786. $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
  787. $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
  788. *
  789. ELSE IF( ITYPE.EQ.8 ) THEN
  790. *
  791. * Symmetric, random eigenvalues
  792. *
  793. CALL SLATMR( N, N, 'S', ISEED, 'S', WORK, 6, ONE, ONE,
  794. $ 'T', 'N', WORK( N+1 ), 1, ONE,
  795. $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
  796. $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
  797. *
  798. ELSE IF( ITYPE.EQ.9 ) THEN
  799. *
  800. * General, random eigenvalues
  801. *
  802. CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
  803. $ 'T', 'N', WORK( N+1 ), 1, ONE,
  804. $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
  805. $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
  806. IF( N.GE.4 ) THEN
  807. CALL SLASET( 'Full', 2, N, ZERO, ZERO, A, LDA )
  808. CALL SLASET( 'Full', N-3, 1, ZERO, ZERO, A( 3, 1 ),
  809. $ LDA )
  810. CALL SLASET( 'Full', N-3, 2, ZERO, ZERO, A( 3, N-1 ),
  811. $ LDA )
  812. CALL SLASET( 'Full', 1, N, ZERO, ZERO, A( N, 1 ),
  813. $ LDA )
  814. END IF
  815. *
  816. ELSE IF( ITYPE.EQ.10 ) THEN
  817. *
  818. * Triangular, random eigenvalues
  819. *
  820. CALL SLATMR( N, N, 'S', ISEED, 'N', WORK, 6, ONE, ONE,
  821. $ 'T', 'N', WORK( N+1 ), 1, ONE,
  822. $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
  823. $ ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
  824. *
  825. ELSE
  826. *
  827. IINFO = 1
  828. END IF
  829. *
  830. IF( IINFO.NE.0 ) THEN
  831. WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
  832. $ IOLDSD
  833. INFO = ABS( IINFO )
  834. RETURN
  835. END IF
  836. *
  837. 90 CONTINUE
  838. *
  839. * Test for minimal and generous workspace
  840. *
  841. DO 130 IWK = 1, 3
  842. IF( IWK.EQ.1 ) THEN
  843. NNWORK = 3*N
  844. ELSE IF( IWK.EQ.2 ) THEN
  845. NNWORK = 6*N + N**2
  846. ELSE
  847. NNWORK = 6*N + 2*N**2
  848. END IF
  849. NNWORK = MAX( NNWORK, 1 )
  850. *
  851. * Test for all balancing options
  852. *
  853. DO 120 IBAL = 1, 4
  854. BALANC = BAL( IBAL )
  855. *
  856. * Perform tests
  857. *
  858. CALL SGET23( .FALSE., BALANC, JTYPE, THRESH, IOLDSD,
  859. $ NOUNIT, N, A, LDA, H, WR, WI, WR1, WI1,
  860. $ VL, LDVL, VR, LDVR, LRE, LDLRE, RCONDV,
  861. $ RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
  862. $ SCALE, SCALE1, RESULT, WORK, NNWORK,
  863. $ IWORK, INFO )
  864. *
  865. * Check for RESULT(j) > THRESH
  866. *
  867. NTEST = 0
  868. NFAIL = 0
  869. DO 100 J = 1, 9
  870. IF( RESULT( J ).GE.ZERO )
  871. $ NTEST = NTEST + 1
  872. IF( RESULT( J ).GE.THRESH )
  873. $ NFAIL = NFAIL + 1
  874. 100 CONTINUE
  875. *
  876. IF( NFAIL.GT.0 )
  877. $ NTESTF = NTESTF + 1
  878. IF( NTESTF.EQ.1 ) THEN
  879. WRITE( NOUNIT, FMT = 9999 )PATH
  880. WRITE( NOUNIT, FMT = 9998 )
  881. WRITE( NOUNIT, FMT = 9997 )
  882. WRITE( NOUNIT, FMT = 9996 )
  883. WRITE( NOUNIT, FMT = 9995 )THRESH
  884. NTESTF = 2
  885. END IF
  886. *
  887. DO 110 J = 1, 9
  888. IF( RESULT( J ).GE.THRESH ) THEN
  889. WRITE( NOUNIT, FMT = 9994 )BALANC, N, IWK,
  890. $ IOLDSD, JTYPE, J, RESULT( J )
  891. END IF
  892. 110 CONTINUE
  893. *
  894. NERRS = NERRS + NFAIL
  895. NTESTT = NTESTT + NTEST
  896. *
  897. 120 CONTINUE
  898. 130 CONTINUE
  899. 140 CONTINUE
  900. 150 CONTINUE
  901. *
  902. 160 CONTINUE
  903. *
  904. * Read in data from file to check accuracy of condition estimation.
  905. * Assume input eigenvalues are sorted lexicographically (increasing
  906. * by real part, then decreasing by imaginary part)
  907. *
  908. JTYPE = 0
  909. 170 CONTINUE
  910. READ( NIUNIT, FMT = *, END = 220 )N
  911. *
  912. * Read input data until N=0
  913. *
  914. IF( N.EQ.0 )
  915. $ GO TO 220
  916. JTYPE = JTYPE + 1
  917. ISEED( 1 ) = JTYPE
  918. DO 180 I = 1, N
  919. READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
  920. 180 CONTINUE
  921. DO 190 I = 1, N
  922. READ( NIUNIT, FMT = * )WR1( I ), WI1( I ), RCDEIN( I ),
  923. $ RCDVIN( I )
  924. 190 CONTINUE
  925. CALL SGET23( .TRUE., 'N', 22, THRESH, ISEED, NOUNIT, N, A, LDA, H,
  926. $ WR, WI, WR1, WI1, VL, LDVL, VR, LDVR, LRE, LDLRE,
  927. $ RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
  928. $ SCALE, SCALE1, RESULT, WORK, 6*N+2*N**2, IWORK,
  929. $ INFO )
  930. *
  931. * Check for RESULT(j) > THRESH
  932. *
  933. NTEST = 0
  934. NFAIL = 0
  935. DO 200 J = 1, 11
  936. IF( RESULT( J ).GE.ZERO )
  937. $ NTEST = NTEST + 1
  938. IF( RESULT( J ).GE.THRESH )
  939. $ NFAIL = NFAIL + 1
  940. 200 CONTINUE
  941. *
  942. IF( NFAIL.GT.0 )
  943. $ NTESTF = NTESTF + 1
  944. IF( NTESTF.EQ.1 ) THEN
  945. WRITE( NOUNIT, FMT = 9999 )PATH
  946. WRITE( NOUNIT, FMT = 9998 )
  947. WRITE( NOUNIT, FMT = 9997 )
  948. WRITE( NOUNIT, FMT = 9996 )
  949. WRITE( NOUNIT, FMT = 9995 )THRESH
  950. NTESTF = 2
  951. END IF
  952. *
  953. DO 210 J = 1, 11
  954. IF( RESULT( J ).GE.THRESH ) THEN
  955. WRITE( NOUNIT, FMT = 9993 )N, JTYPE, J, RESULT( J )
  956. END IF
  957. 210 CONTINUE
  958. *
  959. NERRS = NERRS + NFAIL
  960. NTESTT = NTESTT + NTEST
  961. GO TO 170
  962. 220 CONTINUE
  963. *
  964. * Summary
  965. *
  966. CALL SLASUM( PATH, NOUNIT, NERRS, NTESTT )
  967. *
  968. 9999 FORMAT( / 1X, A3, ' -- Real Eigenvalue-Eigenvector Decomposition',
  969. $ ' Expert Driver', /
  970. $ ' Matrix types (see SDRVVX for details): ' )
  971. *
  972. 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
  973. $ ' ', ' 5=Diagonal: geometr. spaced entries.',
  974. $ / ' 2=Identity matrix. ', ' 6=Diagona',
  975. $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
  976. $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
  977. $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
  978. $ 'mall, evenly spaced.' )
  979. 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
  980. $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
  981. $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
  982. $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
  983. $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
  984. $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
  985. $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
  986. $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
  987. $ ' complx ' )
  988. 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
  989. $ 'with small random entries.', / ' 20=Matrix with large ran',
  990. $ 'dom entries. ', ' 22=Matrix read from input file', / )
  991. 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
  992. $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ',
  993. $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ',
  994. $ / ' 3 = | |VR(i)| - 1 | / ulp ',
  995. $ / ' 4 = | |VL(i)| - 1 | / ulp ',
  996. $ / ' 5 = 0 if W same no matter if VR or VL computed,',
  997. $ ' 1/ulp otherwise', /
  998. $ ' 6 = 0 if VR same no matter what else computed,',
  999. $ ' 1/ulp otherwise', /
  1000. $ ' 7 = 0 if VL same no matter what else computed,',
  1001. $ ' 1/ulp otherwise', /
  1002. $ ' 8 = 0 if RCONDV same no matter what else computed,',
  1003. $ ' 1/ulp otherwise', /
  1004. $ ' 9 = 0 if SCALE, ILO, IHI, ABNRM same no matter what else',
  1005. $ ' computed, 1/ulp otherwise',
  1006. $ / ' 10 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),',
  1007. $ / ' 11 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),' )
  1008. 9994 FORMAT( ' BALANC=''', A1, ''',N=', I4, ',IWK=', I1, ', seed=',
  1009. $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 )
  1010. 9993 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
  1011. $ G10.3 )
  1012. 9992 FORMAT( ' SDRVVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  1013. $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
  1014. *
  1015. RETURN
  1016. *
  1017. * End of SDRVVX
  1018. *
  1019. END