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.

cchkgg.f 45 kB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244
  1. *> \brief \b CCHKGG
  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 CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  12. * TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
  13. * S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
  14. * ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
  15. * RWORK, LLWORK, RESULT, INFO )
  16. *
  17. * .. Scalar Arguments ..
  18. * LOGICAL TSTDIF
  19. * INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
  20. * REAL THRESH, THRSHN
  21. * ..
  22. * .. Array Arguments ..
  23. * LOGICAL DOTYPE( * ), LLWORK( * )
  24. * INTEGER ISEED( 4 ), NN( * )
  25. * REAL RESULT( 15 ), RWORK( * )
  26. * COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
  27. * $ B( LDA, * ), BETA1( * ), BETA3( * ),
  28. * $ EVECTL( LDU, * ), EVECTR( LDU, * ),
  29. * $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
  30. * $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
  31. * $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
  32. * $ WORK( * ), Z( LDU, * )
  33. * ..
  34. *
  35. *
  36. *> \par Purpose:
  37. * =============
  38. *>
  39. *> \verbatim
  40. *>
  41. *> CCHKGG checks the nonsymmetric generalized eigenvalue problem
  42. *> routines.
  43. *> H H H
  44. *> CGGHRD factors A and B as U H V and U T V , where means conjugate
  45. *> transpose, H is hessenberg, T is triangular and U and V are unitary.
  46. *>
  47. *> H H
  48. *> CHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper
  49. *> triangular and Q and Z are unitary. It also computes the generalized
  50. *> eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where
  51. *> alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j)
  52. *> is a root of the generalized eigenvalue problem
  53. *>
  54. *> det( A - w(j) B ) = 0
  55. *>
  56. *> and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
  57. *> problem
  58. *>
  59. *> det( m(j) A - B ) = 0
  60. *>
  61. *> CTGEVC computes the matrix L of left eigenvectors and the matrix R
  62. *> of right eigenvectors for the matrix pair ( S, P ). In the
  63. *> description below, l and r are left and right eigenvectors
  64. *> corresponding to the generalized eigenvalues (alpha,beta).
  65. *>
  66. *> When CCHKGG is called, a number of matrix "sizes" ("n's") and a
  67. *> number of matrix "types" are specified. For each size ("n")
  68. *> and each type of matrix, one matrix will be generated and used
  69. *> to test the nonsymmetric eigenroutines. For each matrix, 13
  70. *> tests will be performed. The first twelve "test ratios" should be
  71. *> small -- O(1). They will be compared with the threshold THRESH:
  72. *>
  73. *> H
  74. *> (1) | A - U H V | / ( |A| n ulp )
  75. *>
  76. *> H
  77. *> (2) | B - U T V | / ( |B| n ulp )
  78. *>
  79. *> H
  80. *> (3) | I - UU | / ( n ulp )
  81. *>
  82. *> H
  83. *> (4) | I - VV | / ( n ulp )
  84. *>
  85. *> H
  86. *> (5) | H - Q S Z | / ( |H| n ulp )
  87. *>
  88. *> H
  89. *> (6) | T - Q P Z | / ( |T| n ulp )
  90. *>
  91. *> H
  92. *> (7) | I - QQ | / ( n ulp )
  93. *>
  94. *> H
  95. *> (8) | I - ZZ | / ( n ulp )
  96. *>
  97. *> (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
  98. *> H
  99. *> | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) )
  100. *>
  101. *> (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of
  102. *> H
  103. *> | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) )
  104. *>
  105. *> where the eigenvectors l' are the result of passing Q to
  106. *> STGEVC and back transforming (JOB='B').
  107. *>
  108. *> (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
  109. *>
  110. *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
  111. *>
  112. *> (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of
  113. *>
  114. *> | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
  115. *>
  116. *> where the eigenvectors r' are the result of passing Z to
  117. *> STGEVC and back transforming (JOB='B').
  118. *>
  119. *> The last three test ratios will usually be small, but there is no
  120. *> mathematical requirement that they be so. They are therefore
  121. *> compared with THRESH only if TSTDIF is .TRUE.
  122. *>
  123. *> (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
  124. *>
  125. *> (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
  126. *>
  127. *> (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
  128. *> |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
  129. *>
  130. *> In addition, the normalization of L and R are checked, and compared
  131. *> with the threshold THRSHN.
  132. *>
  133. *> Test Matrices
  134. *> ---- --------
  135. *>
  136. *> The sizes of the test matrices are specified by an array
  137. *> NN(1:NSIZES); the value of each element NN(j) specifies one size.
  138. *> The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
  139. *> DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
  140. *> Currently, the list of possible types is:
  141. *>
  142. *> (1) ( 0, 0 ) (a pair of zero matrices)
  143. *>
  144. *> (2) ( I, 0 ) (an identity and a zero matrix)
  145. *>
  146. *> (3) ( 0, I ) (an identity and a zero matrix)
  147. *>
  148. *> (4) ( I, I ) (a pair of identity matrices)
  149. *>
  150. *> t t
  151. *> (5) ( J , J ) (a pair of transposed Jordan blocks)
  152. *>
  153. *> t ( I 0 )
  154. *> (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
  155. *> ( 0 I ) ( 0 J )
  156. *> and I is a k x k identity and J a (k+1)x(k+1)
  157. *> Jordan block; k=(N-1)/2
  158. *>
  159. *> (7) ( D, I ) where D is P*D1, P is a random unitary diagonal
  160. *> matrix (i.e., with random magnitude 1 entries
  161. *> on the diagonal), and D1=diag( 0, 1,..., N-1 )
  162. *> (i.e., a diagonal matrix with D1(1,1)=0,
  163. *> D1(2,2)=1, ..., D1(N,N)=N-1.)
  164. *> (8) ( I, D )
  165. *>
  166. *> (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
  167. *>
  168. *> (10) ( small*D, big*I )
  169. *>
  170. *> (11) ( big*I, small*D )
  171. *>
  172. *> (12) ( small*I, big*D )
  173. *>
  174. *> (13) ( big*D, big*I )
  175. *>
  176. *> (14) ( small*D, small*I )
  177. *>
  178. *> (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and
  179. *> D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and
  180. *> P and Q are random unitary diagonal matrices.
  181. *> t t
  182. *> (16) U ( J , J ) V where U and V are random unitary matrices.
  183. *>
  184. *> (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices
  185. *> with random O(1) entries above the diagonal
  186. *> and diagonal entries diag(T1) =
  187. *> P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
  188. *> Q*( 0, N-3, N-4,..., 1, 0, 0 )
  189. *>
  190. *> (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
  191. *> diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
  192. *> s = machine precision.
  193. *>
  194. *> (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
  195. *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
  196. *>
  197. *> N-5
  198. *> (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
  199. *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
  200. *>
  201. *> (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
  202. *> diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
  203. *> where r1,..., r(N-4) are random.
  204. *>
  205. *> (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
  206. *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
  207. *>
  208. *> (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
  209. *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
  210. *>
  211. *> (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
  212. *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
  213. *>
  214. *> (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
  215. *> diag(T2) = ( 0, 1, ..., 1, 0, 0 )
  216. *>
  217. *> (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular
  218. *> matrices.
  219. *> \endverbatim
  220. *
  221. * Arguments:
  222. * ==========
  223. *
  224. *> \param[in] NSIZES
  225. *> \verbatim
  226. *> NSIZES is INTEGER
  227. *> The number of sizes of matrices to use. If it is zero,
  228. *> CCHKGG does nothing. It must be at least zero.
  229. *> \endverbatim
  230. *>
  231. *> \param[in] NN
  232. *> \verbatim
  233. *> NN is INTEGER array, dimension (NSIZES)
  234. *> An array containing the sizes to be used for the matrices.
  235. *> Zero values will be skipped. The values must be at least
  236. *> zero.
  237. *> \endverbatim
  238. *>
  239. *> \param[in] NTYPES
  240. *> \verbatim
  241. *> NTYPES is INTEGER
  242. *> The number of elements in DOTYPE. If it is zero, CCHKGG
  243. *> does nothing. It must be at least zero. If it is MAXTYP+1
  244. *> and NSIZES is 1, then an additional type, MAXTYP+1 is
  245. *> defined, which is to use whatever matrix is in A. This
  246. *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
  247. *> DOTYPE(MAXTYP+1) is .TRUE. .
  248. *> \endverbatim
  249. *>
  250. *> \param[in] DOTYPE
  251. *> \verbatim
  252. *> DOTYPE is LOGICAL array, dimension (NTYPES)
  253. *> If DOTYPE(j) is .TRUE., then for each size in NN a
  254. *> matrix of that size and of type j will be generated.
  255. *> If NTYPES is smaller than the maximum number of types
  256. *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
  257. *> MAXTYP will not be generated. If NTYPES is larger
  258. *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
  259. *> will be ignored.
  260. *> \endverbatim
  261. *>
  262. *> \param[in,out] ISEED
  263. *> \verbatim
  264. *> ISEED is INTEGER array, dimension (4)
  265. *> On entry ISEED specifies the seed of the random number
  266. *> generator. The array elements should be between 0 and 4095;
  267. *> if not they will be reduced mod 4096. Also, ISEED(4) must
  268. *> be odd. The random number generator uses a linear
  269. *> congruential sequence limited to small integers, and so
  270. *> should produce machine independent random numbers. The
  271. *> values of ISEED are changed on exit, and can be used in the
  272. *> next call to CCHKGG to continue the same random number
  273. *> sequence.
  274. *> \endverbatim
  275. *>
  276. *> \param[in] THRESH
  277. *> \verbatim
  278. *> THRESH is REAL
  279. *> A test will count as "failed" if the "error", computed as
  280. *> described above, exceeds THRESH. Note that the error
  281. *> is scaled to be O(1), so THRESH should be a reasonably
  282. *> small multiple of 1, e.g., 10 or 100. In particular,
  283. *> it should not depend on the precision (single vs. double)
  284. *> or the size of the matrix. It must be at least zero.
  285. *> \endverbatim
  286. *>
  287. *> \param[in] TSTDIF
  288. *> \verbatim
  289. *> TSTDIF is LOGICAL
  290. *> Specifies whether test ratios 13-15 will be computed and
  291. *> compared with THRESH.
  292. *> = .FALSE.: Only test ratios 1-12 will be computed and tested.
  293. *> Ratios 13-15 will be set to zero.
  294. *> = .TRUE.: All the test ratios 1-15 will be computed and
  295. *> tested.
  296. *> \endverbatim
  297. *>
  298. *> \param[in] THRSHN
  299. *> \verbatim
  300. *> THRSHN is REAL
  301. *> Threshold for reporting eigenvector normalization error.
  302. *> If the normalization of any eigenvector differs from 1 by
  303. *> more than THRSHN*ulp, then a special error message will be
  304. *> printed. (This is handled separately from the other tests,
  305. *> since only a compiler or programming error should cause an
  306. *> error message, at least if THRSHN is at least 5--10.)
  307. *> \endverbatim
  308. *>
  309. *> \param[in] NOUNIT
  310. *> \verbatim
  311. *> NOUNIT is INTEGER
  312. *> The FORTRAN unit number for printing out error messages
  313. *> (e.g., if a routine returns IINFO not equal to 0.)
  314. *> \endverbatim
  315. *>
  316. *> \param[in,out] A
  317. *> \verbatim
  318. *> A is COMPLEX array, dimension (LDA, max(NN))
  319. *> Used to hold the original A matrix. Used as input only
  320. *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
  321. *> DOTYPE(MAXTYP+1)=.TRUE.
  322. *> \endverbatim
  323. *>
  324. *> \param[in] LDA
  325. *> \verbatim
  326. *> LDA is INTEGER
  327. *> The leading dimension of A, B, H, T, S1, P1, S2, and P2.
  328. *> It must be at least 1 and at least max( NN ).
  329. *> \endverbatim
  330. *>
  331. *> \param[in,out] B
  332. *> \verbatim
  333. *> B is COMPLEX array, dimension (LDA, max(NN))
  334. *> Used to hold the original B matrix. Used as input only
  335. *> if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
  336. *> DOTYPE(MAXTYP+1)=.TRUE.
  337. *> \endverbatim
  338. *>
  339. *> \param[out] H
  340. *> \verbatim
  341. *> H is COMPLEX array, dimension (LDA, max(NN))
  342. *> The upper Hessenberg matrix computed from A by CGGHRD.
  343. *> \endverbatim
  344. *>
  345. *> \param[out] T
  346. *> \verbatim
  347. *> T is COMPLEX array, dimension (LDA, max(NN))
  348. *> The upper triangular matrix computed from B by CGGHRD.
  349. *> \endverbatim
  350. *>
  351. *> \param[out] S1
  352. *> \verbatim
  353. *> S1 is COMPLEX array, dimension (LDA, max(NN))
  354. *> The Schur (upper triangular) matrix computed from H by CHGEQZ
  355. *> when Q and Z are also computed.
  356. *> \endverbatim
  357. *>
  358. *> \param[out] S2
  359. *> \verbatim
  360. *> S2 is COMPLEX array, dimension (LDA, max(NN))
  361. *> The Schur (upper triangular) matrix computed from H by CHGEQZ
  362. *> when Q and Z are not computed.
  363. *> \endverbatim
  364. *>
  365. *> \param[out] P1
  366. *> \verbatim
  367. *> P1 is COMPLEX array, dimension (LDA, max(NN))
  368. *> The upper triangular matrix computed from T by CHGEQZ
  369. *> when Q and Z are also computed.
  370. *> \endverbatim
  371. *>
  372. *> \param[out] P2
  373. *> \verbatim
  374. *> P2 is COMPLEX array, dimension (LDA, max(NN))
  375. *> The upper triangular matrix computed from T by CHGEQZ
  376. *> when Q and Z are not computed.
  377. *> \endverbatim
  378. *>
  379. *> \param[out] U
  380. *> \verbatim
  381. *> U is COMPLEX array, dimension (LDU, max(NN))
  382. *> The (left) unitary matrix computed by CGGHRD.
  383. *> \endverbatim
  384. *>
  385. *> \param[in] LDU
  386. *> \verbatim
  387. *> LDU is INTEGER
  388. *> The leading dimension of U, V, Q, Z, EVECTL, and EVECTR. It
  389. *> must be at least 1 and at least max( NN ).
  390. *> \endverbatim
  391. *>
  392. *> \param[out] V
  393. *> \verbatim
  394. *> V is COMPLEX array, dimension (LDU, max(NN))
  395. *> The (right) unitary matrix computed by CGGHRD.
  396. *> \endverbatim
  397. *>
  398. *> \param[out] Q
  399. *> \verbatim
  400. *> Q is COMPLEX array, dimension (LDU, max(NN))
  401. *> The (left) unitary matrix computed by CHGEQZ.
  402. *> \endverbatim
  403. *>
  404. *> \param[out] Z
  405. *> \verbatim
  406. *> Z is COMPLEX array, dimension (LDU, max(NN))
  407. *> The (left) unitary matrix computed by CHGEQZ.
  408. *> \endverbatim
  409. *>
  410. *> \param[out] ALPHA1
  411. *> \verbatim
  412. *> ALPHA1 is COMPLEX array, dimension (max(NN))
  413. *> \endverbatim
  414. *>
  415. *> \param[out] BETA1
  416. *> \verbatim
  417. *> BETA1 is COMPLEX array, dimension (max(NN))
  418. *> The generalized eigenvalues of (A,B) computed by CHGEQZ
  419. *> when Q, Z, and the full Schur matrices are computed.
  420. *> \endverbatim
  421. *>
  422. *> \param[out] ALPHA3
  423. *> \verbatim
  424. *> ALPHA3 is COMPLEX array, dimension (max(NN))
  425. *> \endverbatim
  426. *>
  427. *> \param[out] BETA3
  428. *> \verbatim
  429. *> BETA3 is COMPLEX array, dimension (max(NN))
  430. *> The generalized eigenvalues of (A,B) computed by CHGEQZ
  431. *> when neither Q, Z, nor the Schur matrices are computed.
  432. *> \endverbatim
  433. *>
  434. *> \param[out] EVECTL
  435. *> \verbatim
  436. *> EVECTL is COMPLEX array, dimension (LDU, max(NN))
  437. *> The (lower triangular) left eigenvector matrix for the
  438. *> matrices in S1 and P1.
  439. *> \endverbatim
  440. *>
  441. *> \param[out] EVECTR
  442. *> \verbatim
  443. *> EVECTR is COMPLEX array, dimension (LDU, max(NN))
  444. *> The (upper triangular) right eigenvector matrix for the
  445. *> matrices in S1 and P1.
  446. *> \endverbatim
  447. *>
  448. *> \param[out] WORK
  449. *> \verbatim
  450. *> WORK is COMPLEX array, dimension (LWORK)
  451. *> \endverbatim
  452. *>
  453. *> \param[in] LWORK
  454. *> \verbatim
  455. *> LWORK is INTEGER
  456. *> The number of entries in WORK. This must be at least
  457. *> max( 4*N, 2 * N**2, 1 ), for all N=NN(j).
  458. *> \endverbatim
  459. *>
  460. *> \param[out] RWORK
  461. *> \verbatim
  462. *> RWORK is REAL array, dimension (2*max(NN))
  463. *> \endverbatim
  464. *>
  465. *> \param[out] LLWORK
  466. *> \verbatim
  467. *> LLWORK is LOGICAL array, dimension (max(NN))
  468. *> \endverbatim
  469. *>
  470. *> \param[out] RESULT
  471. *> \verbatim
  472. *> RESULT is REAL array, dimension (15)
  473. *> The values computed by the tests described above.
  474. *> The values are currently limited to 1/ulp, to avoid
  475. *> overflow.
  476. *> \endverbatim
  477. *>
  478. *> \param[out] INFO
  479. *> \verbatim
  480. *> INFO is INTEGER
  481. *> = 0: successful exit.
  482. *> < 0: if INFO = -i, the i-th argument had an illegal value.
  483. *> > 0: A routine returned an error code. INFO is the
  484. *> absolute value of the INFO value returned.
  485. *> \endverbatim
  486. *
  487. * Authors:
  488. * ========
  489. *
  490. *> \author Univ. of Tennessee
  491. *> \author Univ. of California Berkeley
  492. *> \author Univ. of Colorado Denver
  493. *> \author NAG Ltd.
  494. *
  495. *> \date June 2016
  496. *
  497. *> \ingroup complex_eig
  498. *
  499. * =====================================================================
  500. SUBROUTINE CCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
  501. $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
  502. $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
  503. $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
  504. $ RWORK, LLWORK, RESULT, INFO )
  505. *
  506. * -- LAPACK test routine (version 3.7.0) --
  507. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  508. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  509. * June 2016
  510. *
  511. * .. Scalar Arguments ..
  512. LOGICAL TSTDIF
  513. INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
  514. REAL THRESH, THRSHN
  515. * ..
  516. * .. Array Arguments ..
  517. LOGICAL DOTYPE( * ), LLWORK( * )
  518. INTEGER ISEED( 4 ), NN( * )
  519. REAL RESULT( 15 ), RWORK( * )
  520. COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
  521. $ B( LDA, * ), BETA1( * ), BETA3( * ),
  522. $ EVECTL( LDU, * ), EVECTR( LDU, * ),
  523. $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
  524. $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
  525. $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
  526. $ WORK( * ), Z( LDU, * )
  527. * ..
  528. *
  529. * =====================================================================
  530. *
  531. * .. Parameters ..
  532. REAL ZERO, ONE
  533. PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 )
  534. COMPLEX CZERO, CONE
  535. PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ),
  536. $ CONE = ( 1.0E+0, 0.0E+0 ) )
  537. INTEGER MAXTYP
  538. PARAMETER ( MAXTYP = 26 )
  539. * ..
  540. * .. Local Scalars ..
  541. LOGICAL BADNN
  542. INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
  543. $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
  544. $ NTEST, NTESTT
  545. REAL ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
  546. $ ULP, ULPINV
  547. COMPLEX CTEMP
  548. * ..
  549. * .. Local Arrays ..
  550. LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
  551. INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
  552. $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
  553. $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
  554. $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
  555. $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
  556. REAL DUMMA( 4 ), RMAGN( 0: 3 )
  557. COMPLEX CDUMMA( 4 )
  558. * ..
  559. * .. External Functions ..
  560. REAL CLANGE, SLAMCH
  561. COMPLEX CLARND
  562. EXTERNAL CLANGE, SLAMCH, CLARND
  563. * ..
  564. * .. External Subroutines ..
  565. EXTERNAL CGEQR2, CGET51, CGET52, CGGHRD, CHGEQZ, CLACPY,
  566. $ CLARFG, CLASET, CLATM4, CTGEVC, CUNM2R, SLABAD,
  567. $ SLASUM, XERBLA
  568. * ..
  569. * .. Intrinsic Functions ..
  570. INTRINSIC ABS, CONJG, MAX, MIN, REAL, SIGN
  571. * ..
  572. * .. Data statements ..
  573. DATA KCLASS / 15*1, 10*2, 1*3 /
  574. DATA KZ1 / 0, 1, 2, 1, 3, 3 /
  575. DATA KZ2 / 0, 0, 1, 2, 1, 1 /
  576. DATA KADD / 0, 0, 0, 0, 3, 2 /
  577. DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
  578. $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
  579. DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
  580. $ 1, 1, -4, 2, -4, 8*8, 0 /
  581. DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
  582. $ 4*5, 4*3, 1 /
  583. DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
  584. $ 4*6, 4*4, 1 /
  585. DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
  586. $ 2, 1 /
  587. DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
  588. $ 2, 1 /
  589. DATA KTRIAN / 16*0, 10*1 /
  590. DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
  591. $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
  592. $ 3*.FALSE., 5*.TRUE., .FALSE. /
  593. DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
  594. $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
  595. $ 9*.FALSE. /
  596. * ..
  597. * .. Executable Statements ..
  598. *
  599. * Check for errors
  600. *
  601. INFO = 0
  602. *
  603. BADNN = .FALSE.
  604. NMAX = 1
  605. DO 10 J = 1, NSIZES
  606. NMAX = MAX( NMAX, NN( J ) )
  607. IF( NN( J ).LT.0 )
  608. $ BADNN = .TRUE.
  609. 10 CONTINUE
  610. *
  611. LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 )
  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.LE.1 .OR. LDA.LT.NMAX ) THEN
  624. INFO = -10
  625. ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
  626. INFO = -19
  627. ELSE IF( LWKOPT.GT.LWORK ) THEN
  628. INFO = -30
  629. END IF
  630. *
  631. IF( INFO.NE.0 ) THEN
  632. CALL XERBLA( 'CCHKGG', -INFO )
  633. RETURN
  634. END IF
  635. *
  636. * Quick return if possible
  637. *
  638. IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
  639. $ RETURN
  640. *
  641. SAFMIN = SLAMCH( 'Safe minimum' )
  642. ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
  643. SAFMIN = SAFMIN / ULP
  644. SAFMAX = ONE / SAFMIN
  645. CALL SLABAD( SAFMIN, SAFMAX )
  646. ULPINV = ONE / ULP
  647. *
  648. * The values RMAGN(2:3) depend on N, see below.
  649. *
  650. RMAGN( 0 ) = ZERO
  651. RMAGN( 1 ) = ONE
  652. *
  653. * Loop over sizes, types
  654. *
  655. NTESTT = 0
  656. NERRS = 0
  657. NMATS = 0
  658. *
  659. DO 240 JSIZE = 1, NSIZES
  660. N = NN( JSIZE )
  661. N1 = MAX( 1, N )
  662. RMAGN( 2 ) = SAFMAX*ULP / REAL( N1 )
  663. RMAGN( 3 ) = SAFMIN*ULPINV*N1
  664. *
  665. IF( NSIZES.NE.1 ) THEN
  666. MTYPES = MIN( MAXTYP, NTYPES )
  667. ELSE
  668. MTYPES = MIN( MAXTYP+1, NTYPES )
  669. END IF
  670. *
  671. DO 230 JTYPE = 1, MTYPES
  672. IF( .NOT.DOTYPE( JTYPE ) )
  673. $ GO TO 230
  674. NMATS = NMATS + 1
  675. NTEST = 0
  676. *
  677. * Save ISEED in case of an error.
  678. *
  679. DO 20 J = 1, 4
  680. IOLDSD( J ) = ISEED( J )
  681. 20 CONTINUE
  682. *
  683. * Initialize RESULT
  684. *
  685. DO 30 J = 1, 15
  686. RESULT( J ) = ZERO
  687. 30 CONTINUE
  688. *
  689. * Compute A and B
  690. *
  691. * Description of control parameters:
  692. *
  693. * KCLASS: =1 means w/o rotation, =2 means w/ rotation,
  694. * =3 means random.
  695. * KATYPE: the "type" to be passed to CLATM4 for computing A.
  696. * KAZERO: the pattern of zeros on the diagonal for A:
  697. * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
  698. * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
  699. * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
  700. * non-zero entries.)
  701. * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
  702. * =2: large, =3: small.
  703. * LASIGN: .TRUE. if the diagonal elements of A are to be
  704. * multiplied by a random magnitude 1 number.
  705. * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
  706. * KTRIAN: =0: don't fill in the upper triangle, =1: do.
  707. * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
  708. * RMAGN: used to implement KAMAGN and KBMAGN.
  709. *
  710. IF( MTYPES.GT.MAXTYP )
  711. $ GO TO 110
  712. IINFO = 0
  713. IF( KCLASS( JTYPE ).LT.3 ) THEN
  714. *
  715. * Generate A (w/o rotation)
  716. *
  717. IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
  718. IN = 2*( ( N-1 ) / 2 ) + 1
  719. IF( IN.NE.N )
  720. $ CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
  721. ELSE
  722. IN = N
  723. END IF
  724. CALL CLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
  725. $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
  726. $ RMAGN( KAMAGN( JTYPE ) ), ULP,
  727. $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4,
  728. $ ISEED, A, LDA )
  729. IADD = KADD( KAZERO( JTYPE ) )
  730. IF( IADD.GT.0 .AND. IADD.LE.N )
  731. $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
  732. *
  733. * Generate B (w/o rotation)
  734. *
  735. IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
  736. IN = 2*( ( N-1 ) / 2 ) + 1
  737. IF( IN.NE.N )
  738. $ CALL CLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
  739. ELSE
  740. IN = N
  741. END IF
  742. CALL CLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
  743. $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
  744. $ RMAGN( KBMAGN( JTYPE ) ), ONE,
  745. $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4,
  746. $ ISEED, B, LDA )
  747. IADD = KADD( KBZERO( JTYPE ) )
  748. IF( IADD.NE.0 )
  749. $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
  750. *
  751. IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
  752. *
  753. * Include rotations
  754. *
  755. * Generate U, V as Householder transformations times a
  756. * diagonal matrix. (Note that CLARFG makes U(j,j) and
  757. * V(j,j) real.)
  758. *
  759. DO 50 JC = 1, N - 1
  760. DO 40 JR = JC, N
  761. U( JR, JC ) = CLARND( 3, ISEED )
  762. V( JR, JC ) = CLARND( 3, ISEED )
  763. 40 CONTINUE
  764. CALL CLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
  765. $ WORK( JC ) )
  766. WORK( 2*N+JC ) = SIGN( ONE, REAL( U( JC, JC ) ) )
  767. U( JC, JC ) = CONE
  768. CALL CLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
  769. $ WORK( N+JC ) )
  770. WORK( 3*N+JC ) = SIGN( ONE, REAL( V( JC, JC ) ) )
  771. V( JC, JC ) = CONE
  772. 50 CONTINUE
  773. CTEMP = CLARND( 3, ISEED )
  774. U( N, N ) = CONE
  775. WORK( N ) = CZERO
  776. WORK( 3*N ) = CTEMP / ABS( CTEMP )
  777. CTEMP = CLARND( 3, ISEED )
  778. V( N, N ) = CONE
  779. WORK( 2*N ) = CZERO
  780. WORK( 4*N ) = CTEMP / ABS( CTEMP )
  781. *
  782. * Apply the diagonal matrices
  783. *
  784. DO 70 JC = 1, N
  785. DO 60 JR = 1, N
  786. A( JR, JC ) = WORK( 2*N+JR )*
  787. $ CONJG( WORK( 3*N+JC ) )*
  788. $ A( JR, JC )
  789. B( JR, JC ) = WORK( 2*N+JR )*
  790. $ CONJG( WORK( 3*N+JC ) )*
  791. $ B( JR, JC )
  792. 60 CONTINUE
  793. 70 CONTINUE
  794. CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A,
  795. $ LDA, WORK( 2*N+1 ), IINFO )
  796. IF( IINFO.NE.0 )
  797. $ GO TO 100
  798. CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
  799. $ A, LDA, WORK( 2*N+1 ), IINFO )
  800. IF( IINFO.NE.0 )
  801. $ GO TO 100
  802. CALL CUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B,
  803. $ LDA, WORK( 2*N+1 ), IINFO )
  804. IF( IINFO.NE.0 )
  805. $ GO TO 100
  806. CALL CUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
  807. $ B, LDA, WORK( 2*N+1 ), IINFO )
  808. IF( IINFO.NE.0 )
  809. $ GO TO 100
  810. END IF
  811. ELSE
  812. *
  813. * Random matrices
  814. *
  815. DO 90 JC = 1, N
  816. DO 80 JR = 1, N
  817. A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
  818. $ CLARND( 4, ISEED )
  819. B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
  820. $ CLARND( 4, ISEED )
  821. 80 CONTINUE
  822. 90 CONTINUE
  823. END IF
  824. *
  825. ANORM = CLANGE( '1', N, N, A, LDA, RWORK )
  826. BNORM = CLANGE( '1', N, N, B, LDA, RWORK )
  827. *
  828. 100 CONTINUE
  829. *
  830. IF( IINFO.NE.0 ) THEN
  831. WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
  832. $ IOLDSD
  833. INFO = ABS( IINFO )
  834. RETURN
  835. END IF
  836. *
  837. 110 CONTINUE
  838. *
  839. * Call CGEQR2, CUNM2R, and CGGHRD to compute H, T, U, and V
  840. *
  841. CALL CLACPY( ' ', N, N, A, LDA, H, LDA )
  842. CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
  843. NTEST = 1
  844. RESULT( 1 ) = ULPINV
  845. *
  846. CALL CGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
  847. IF( IINFO.NE.0 ) THEN
  848. WRITE( NOUNIT, FMT = 9999 )'CGEQR2', IINFO, N, JTYPE,
  849. $ IOLDSD
  850. INFO = ABS( IINFO )
  851. GO TO 210
  852. END IF
  853. *
  854. CALL CUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA,
  855. $ WORK( N+1 ), IINFO )
  856. IF( IINFO.NE.0 ) THEN
  857. WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE,
  858. $ IOLDSD
  859. INFO = ABS( IINFO )
  860. GO TO 210
  861. END IF
  862. *
  863. CALL CLASET( 'Full', N, N, CZERO, CONE, U, LDU )
  864. CALL CUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU,
  865. $ WORK( N+1 ), IINFO )
  866. IF( IINFO.NE.0 ) THEN
  867. WRITE( NOUNIT, FMT = 9999 )'CUNM2R', IINFO, N, JTYPE,
  868. $ IOLDSD
  869. INFO = ABS( IINFO )
  870. GO TO 210
  871. END IF
  872. *
  873. CALL CGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V,
  874. $ LDU, IINFO )
  875. IF( IINFO.NE.0 ) THEN
  876. WRITE( NOUNIT, FMT = 9999 )'CGGHRD', IINFO, N, JTYPE,
  877. $ IOLDSD
  878. INFO = ABS( IINFO )
  879. GO TO 210
  880. END IF
  881. NTEST = 4
  882. *
  883. * Do tests 1--4
  884. *
  885. CALL CGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
  886. $ RWORK, RESULT( 1 ) )
  887. CALL CGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
  888. $ RWORK, RESULT( 2 ) )
  889. CALL CGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
  890. $ RWORK, RESULT( 3 ) )
  891. CALL CGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
  892. $ RWORK, RESULT( 4 ) )
  893. *
  894. * Call CHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
  895. *
  896. * Compute T1 and UZ
  897. *
  898. * Eigenvalues only
  899. *
  900. CALL CLACPY( ' ', N, N, H, LDA, S2, LDA )
  901. CALL CLACPY( ' ', N, N, T, LDA, P2, LDA )
  902. NTEST = 5
  903. RESULT( 5 ) = ULPINV
  904. *
  905. CALL CHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
  906. $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
  907. $ RWORK, IINFO )
  908. IF( IINFO.NE.0 ) THEN
  909. WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(E)', IINFO, N, JTYPE,
  910. $ IOLDSD
  911. INFO = ABS( IINFO )
  912. GO TO 210
  913. END IF
  914. *
  915. * Eigenvalues and Full Schur Form
  916. *
  917. CALL CLACPY( ' ', N, N, H, LDA, S2, LDA )
  918. CALL CLACPY( ' ', N, N, T, LDA, P2, LDA )
  919. *
  920. CALL CHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
  921. $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
  922. $ RWORK, IINFO )
  923. IF( IINFO.NE.0 ) THEN
  924. WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(S)', IINFO, N, JTYPE,
  925. $ IOLDSD
  926. INFO = ABS( IINFO )
  927. GO TO 210
  928. END IF
  929. *
  930. * Eigenvalues, Schur Form, and Schur Vectors
  931. *
  932. CALL CLACPY( ' ', N, N, H, LDA, S1, LDA )
  933. CALL CLACPY( ' ', N, N, T, LDA, P1, LDA )
  934. *
  935. CALL CHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
  936. $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
  937. $ RWORK, IINFO )
  938. IF( IINFO.NE.0 ) THEN
  939. WRITE( NOUNIT, FMT = 9999 )'CHGEQZ(V)', IINFO, N, JTYPE,
  940. $ IOLDSD
  941. INFO = ABS( IINFO )
  942. GO TO 210
  943. END IF
  944. *
  945. NTEST = 8
  946. *
  947. * Do Tests 5--8
  948. *
  949. CALL CGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
  950. $ RWORK, RESULT( 5 ) )
  951. CALL CGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
  952. $ RWORK, RESULT( 6 ) )
  953. CALL CGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
  954. $ RWORK, RESULT( 7 ) )
  955. CALL CGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
  956. $ RWORK, RESULT( 8 ) )
  957. *
  958. * Compute the Left and Right Eigenvectors of (S1,P1)
  959. *
  960. * 9: Compute the left eigenvector Matrix without
  961. * back transforming:
  962. *
  963. NTEST = 9
  964. RESULT( 9 ) = ULPINV
  965. *
  966. * To test "SELECT" option, compute half of the eigenvectors
  967. * in one call, and half in another
  968. *
  969. I1 = N / 2
  970. DO 120 J = 1, I1
  971. LLWORK( J ) = .TRUE.
  972. 120 CONTINUE
  973. DO 130 J = I1 + 1, N
  974. LLWORK( J ) = .FALSE.
  975. 130 CONTINUE
  976. *
  977. CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
  978. $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
  979. IF( IINFO.NE.0 ) THEN
  980. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S1)', IINFO, N,
  981. $ JTYPE, IOLDSD
  982. INFO = ABS( IINFO )
  983. GO TO 210
  984. END IF
  985. *
  986. I1 = IN
  987. DO 140 J = 1, I1
  988. LLWORK( J ) = .FALSE.
  989. 140 CONTINUE
  990. DO 150 J = I1 + 1, N
  991. LLWORK( J ) = .TRUE.
  992. 150 CONTINUE
  993. *
  994. CALL CTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
  995. $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
  996. $ WORK, RWORK, IINFO )
  997. IF( IINFO.NE.0 ) THEN
  998. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,S2)', IINFO, N,
  999. $ JTYPE, IOLDSD
  1000. INFO = ABS( IINFO )
  1001. GO TO 210
  1002. END IF
  1003. *
  1004. CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
  1005. $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
  1006. RESULT( 9 ) = DUMMA( 1 )
  1007. IF( DUMMA( 2 ).GT.THRSHN ) THEN
  1008. WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=S)',
  1009. $ DUMMA( 2 ), N, JTYPE, IOLDSD
  1010. END IF
  1011. *
  1012. * 10: Compute the left eigenvector Matrix with
  1013. * back transforming:
  1014. *
  1015. NTEST = 10
  1016. RESULT( 10 ) = ULPINV
  1017. CALL CLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
  1018. CALL CTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
  1019. $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
  1020. IF( IINFO.NE.0 ) THEN
  1021. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(L,B)', IINFO, N,
  1022. $ JTYPE, IOLDSD
  1023. INFO = ABS( IINFO )
  1024. GO TO 210
  1025. END IF
  1026. *
  1027. CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
  1028. $ BETA1, WORK, RWORK, DUMMA( 1 ) )
  1029. RESULT( 10 ) = DUMMA( 1 )
  1030. IF( DUMMA( 2 ).GT.THRSHN ) THEN
  1031. WRITE( NOUNIT, FMT = 9998 )'Left', 'CTGEVC(HOWMNY=B)',
  1032. $ DUMMA( 2 ), N, JTYPE, IOLDSD
  1033. END IF
  1034. *
  1035. * 11: Compute the right eigenvector Matrix without
  1036. * back transforming:
  1037. *
  1038. NTEST = 11
  1039. RESULT( 11 ) = ULPINV
  1040. *
  1041. * To test "SELECT" option, compute half of the eigenvectors
  1042. * in one call, and half in another
  1043. *
  1044. I1 = N / 2
  1045. DO 160 J = 1, I1
  1046. LLWORK( J ) = .TRUE.
  1047. 160 CONTINUE
  1048. DO 170 J = I1 + 1, N
  1049. LLWORK( J ) = .FALSE.
  1050. 170 CONTINUE
  1051. *
  1052. CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
  1053. $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
  1054. IF( IINFO.NE.0 ) THEN
  1055. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S1)', IINFO, N,
  1056. $ JTYPE, IOLDSD
  1057. INFO = ABS( IINFO )
  1058. GO TO 210
  1059. END IF
  1060. *
  1061. I1 = IN
  1062. DO 180 J = 1, I1
  1063. LLWORK( J ) = .FALSE.
  1064. 180 CONTINUE
  1065. DO 190 J = I1 + 1, N
  1066. LLWORK( J ) = .TRUE.
  1067. 190 CONTINUE
  1068. *
  1069. CALL CTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
  1070. $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
  1071. $ RWORK, IINFO )
  1072. IF( IINFO.NE.0 ) THEN
  1073. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,S2)', IINFO, N,
  1074. $ JTYPE, IOLDSD
  1075. INFO = ABS( IINFO )
  1076. GO TO 210
  1077. END IF
  1078. *
  1079. CALL CGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
  1080. $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
  1081. RESULT( 11 ) = DUMMA( 1 )
  1082. IF( DUMMA( 2 ).GT.THRESH ) THEN
  1083. WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=S)',
  1084. $ DUMMA( 2 ), N, JTYPE, IOLDSD
  1085. END IF
  1086. *
  1087. * 12: Compute the right eigenvector Matrix with
  1088. * back transforming:
  1089. *
  1090. NTEST = 12
  1091. RESULT( 12 ) = ULPINV
  1092. CALL CLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
  1093. CALL CTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
  1094. $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
  1095. IF( IINFO.NE.0 ) THEN
  1096. WRITE( NOUNIT, FMT = 9999 )'CTGEVC(R,B)', IINFO, N,
  1097. $ JTYPE, IOLDSD
  1098. INFO = ABS( IINFO )
  1099. GO TO 210
  1100. END IF
  1101. *
  1102. CALL CGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
  1103. $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
  1104. RESULT( 12 ) = DUMMA( 1 )
  1105. IF( DUMMA( 2 ).GT.THRESH ) THEN
  1106. WRITE( NOUNIT, FMT = 9998 )'Right', 'CTGEVC(HOWMNY=B)',
  1107. $ DUMMA( 2 ), N, JTYPE, IOLDSD
  1108. END IF
  1109. *
  1110. * Tests 13--15 are done only on request
  1111. *
  1112. IF( TSTDIF ) THEN
  1113. *
  1114. * Do Tests 13--14
  1115. *
  1116. CALL CGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
  1117. $ WORK, RWORK, RESULT( 13 ) )
  1118. CALL CGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
  1119. $ WORK, RWORK, RESULT( 14 ) )
  1120. *
  1121. * Do Test 15
  1122. *
  1123. TEMP1 = ZERO
  1124. TEMP2 = ZERO
  1125. DO 200 J = 1, N
  1126. TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) )
  1127. TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
  1128. 200 CONTINUE
  1129. *
  1130. TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
  1131. TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
  1132. RESULT( 15 ) = MAX( TEMP1, TEMP2 )
  1133. NTEST = 15
  1134. ELSE
  1135. RESULT( 13 ) = ZERO
  1136. RESULT( 14 ) = ZERO
  1137. RESULT( 15 ) = ZERO
  1138. NTEST = 12
  1139. END IF
  1140. *
  1141. * End of Loop -- Check for RESULT(j) > THRESH
  1142. *
  1143. 210 CONTINUE
  1144. *
  1145. NTESTT = NTESTT + NTEST
  1146. *
  1147. * Print out tests which fail.
  1148. *
  1149. DO 220 JR = 1, NTEST
  1150. IF( RESULT( JR ).GE.THRESH ) THEN
  1151. *
  1152. * If this is the first test to fail,
  1153. * print a header to the data file.
  1154. *
  1155. IF( NERRS.EQ.0 ) THEN
  1156. WRITE( NOUNIT, FMT = 9997 )'CGG'
  1157. *
  1158. * Matrix types
  1159. *
  1160. WRITE( NOUNIT, FMT = 9996 )
  1161. WRITE( NOUNIT, FMT = 9995 )
  1162. WRITE( NOUNIT, FMT = 9994 )'Unitary'
  1163. *
  1164. * Tests performed
  1165. *
  1166. WRITE( NOUNIT, FMT = 9993 )'unitary', '*',
  1167. $ 'conjugate transpose', ( '*', J = 1, 10 )
  1168. *
  1169. END IF
  1170. NERRS = NERRS + 1
  1171. IF( RESULT( JR ).LT.10000.0 ) THEN
  1172. WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
  1173. $ RESULT( JR )
  1174. ELSE
  1175. WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
  1176. $ RESULT( JR )
  1177. END IF
  1178. END IF
  1179. 220 CONTINUE
  1180. *
  1181. 230 CONTINUE
  1182. 240 CONTINUE
  1183. *
  1184. * Summary
  1185. *
  1186. CALL SLASUM( 'CGG', NOUNIT, NERRS, NTESTT )
  1187. RETURN
  1188. *
  1189. 9999 FORMAT( ' CCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
  1190. $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
  1191. *
  1192. 9998 FORMAT( ' CCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
  1193. $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
  1194. $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
  1195. $ ')' )
  1196. *
  1197. 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' )
  1198. *
  1199. 9996 FORMAT( ' Matrix types (see CCHKGG for details): ' )
  1200. *
  1201. 9995 FORMAT( ' Special Matrices:', 23X,
  1202. $ '(J''=transposed Jordan block)',
  1203. $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
  1204. $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
  1205. $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
  1206. $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
  1207. $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
  1208. $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
  1209. 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
  1210. $ / ' 16=Transposed Jordan Blocks 19=geometric ',
  1211. $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
  1212. $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
  1213. $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
  1214. $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
  1215. $ '23=(small,large) 24=(small,small) 25=(large,large)',
  1216. $ / ' 26=random O(1) matrices.' )
  1217. *
  1218. 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ',
  1219. $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
  1220. $ ', l and r are the', / 20X,
  1221. $ 'appropriate left and right eigenvectors, resp., a is',
  1222. $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
  1223. $ / ' 1 = | A - U H V', A,
  1224. $ ' | / ( |A| n ulp ) 2 = | B - U T V', A,
  1225. $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
  1226. $ ' | / ( n ulp ) 4 = | I - VV', A,
  1227. $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
  1228. $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
  1229. $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
  1230. $ ' | / ( n ulp ) 8 = | I - ZZ', A,
  1231. $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
  1232. $ ' l | / const. 10 = max | ( b H - a T )', A,
  1233. $ ' l | / const.', /
  1234. $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
  1235. $ ' - a T ) r | / const.', / 1X )
  1236. *
  1237. 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
  1238. $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
  1239. 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
  1240. $ 4( I4, ',' ), ' result ', I2, ' is', 1P, E10.3 )
  1241. *
  1242. * End of CCHKGG
  1243. *
  1244. END