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.

clattp.f 25 kB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822
  1. *> \brief \b CLATTP
  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 CLATTP( IMAT, UPLO, TRANS, DIAG, ISEED, N, AP, B, WORK,
  12. * RWORK, INFO )
  13. *
  14. * .. Scalar Arguments ..
  15. * CHARACTER DIAG, TRANS, UPLO
  16. * INTEGER IMAT, INFO, N
  17. * ..
  18. * .. Array Arguments ..
  19. * INTEGER ISEED( 4 )
  20. * REAL RWORK( * )
  21. * COMPLEX AP( * ), B( * ), WORK( * )
  22. * ..
  23. *
  24. *
  25. *> \par Purpose:
  26. * =============
  27. *>
  28. *> \verbatim
  29. *>
  30. *> CLATTP generates a triangular test matrix in packed storage.
  31. *> IMAT and UPLO uniquely specify the properties of the test matrix,
  32. *> which is returned in the array AP.
  33. *> \endverbatim
  34. *
  35. * Arguments:
  36. * ==========
  37. *
  38. *> \param[in] IMAT
  39. *> \verbatim
  40. *> IMAT is INTEGER
  41. *> An integer key describing which matrix to generate for this
  42. *> path.
  43. *> \endverbatim
  44. *>
  45. *> \param[in] UPLO
  46. *> \verbatim
  47. *> UPLO is CHARACTER*1
  48. *> Specifies whether the matrix A will be upper or lower
  49. *> triangular.
  50. *> = 'U': Upper triangular
  51. *> = 'L': Lower triangular
  52. *> \endverbatim
  53. *>
  54. *> \param[in] TRANS
  55. *> \verbatim
  56. *> TRANS is CHARACTER*1
  57. *> Specifies whether the matrix or its transpose will be used.
  58. *> = 'N': No transpose
  59. *> = 'T': Transpose
  60. *> = 'C': Conjugate transpose
  61. *> \endverbatim
  62. *>
  63. *> \param[out] DIAG
  64. *> \verbatim
  65. *> DIAG is CHARACTER*1
  66. *> Specifies whether or not the matrix A is unit triangular.
  67. *> = 'N': Non-unit triangular
  68. *> = 'U': Unit triangular
  69. *> \endverbatim
  70. *>
  71. *> \param[in,out] ISEED
  72. *> \verbatim
  73. *> ISEED is INTEGER array, dimension (4)
  74. *> The seed vector for the random number generator (used in
  75. *> CLATMS). Modified on exit.
  76. *> \endverbatim
  77. *>
  78. *> \param[in] N
  79. *> \verbatim
  80. *> N is INTEGER
  81. *> The order of the matrix to be generated.
  82. *> \endverbatim
  83. *>
  84. *> \param[out] AP
  85. *> \verbatim
  86. *> AP is COMPLEX array, dimension (N*(N+1)/2)
  87. *> The upper or lower triangular matrix A, packed columnwise in
  88. *> a linear array. The j-th column of A is stored in the array
  89. *> AP as follows:
  90. *> if UPLO = 'U', AP((j-1)*j/2 + i) = A(i,j) for 1<=i<=j;
  91. *> if UPLO = 'L',
  92. *> AP((j-1)*(n-j) + j*(j+1)/2 + i-j) = A(i,j) for j<=i<=n.
  93. *> \endverbatim
  94. *>
  95. *> \param[out] B
  96. *> \verbatim
  97. *> B is COMPLEX array, dimension (N)
  98. *> The right hand side vector, if IMAT > 10.
  99. *> \endverbatim
  100. *>
  101. *> \param[out] WORK
  102. *> \verbatim
  103. *> WORK is COMPLEX array, dimension (2*N)
  104. *> \endverbatim
  105. *>
  106. *> \param[out] RWORK
  107. *> \verbatim
  108. *> RWORK is REAL array, dimension (N)
  109. *> \endverbatim
  110. *>
  111. *> \param[out] INFO
  112. *> \verbatim
  113. *> INFO is INTEGER
  114. *> = 0: successful exit
  115. *> < 0: if INFO = -i, the i-th argument had an illegal value
  116. *> \endverbatim
  117. *
  118. * Authors:
  119. * ========
  120. *
  121. *> \author Univ. of Tennessee
  122. *> \author Univ. of California Berkeley
  123. *> \author Univ. of Colorado Denver
  124. *> \author NAG Ltd.
  125. *
  126. *> \ingroup complex_lin
  127. *
  128. * =====================================================================
  129. SUBROUTINE CLATTP( IMAT, UPLO, TRANS, DIAG, ISEED, N, AP, B, WORK,
  130. $ RWORK, INFO )
  131. *
  132. * -- LAPACK test routine --
  133. * -- LAPACK is a software package provided by Univ. of Tennessee, --
  134. * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  135. *
  136. * .. Scalar Arguments ..
  137. CHARACTER DIAG, TRANS, UPLO
  138. INTEGER IMAT, INFO, N
  139. * ..
  140. * .. Array Arguments ..
  141. INTEGER ISEED( 4 )
  142. REAL RWORK( * )
  143. COMPLEX AP( * ), B( * ), WORK( * )
  144. * ..
  145. *
  146. * =====================================================================
  147. *
  148. * .. Parameters ..
  149. REAL ONE, TWO, ZERO
  150. PARAMETER ( ONE = 1.0E+0, TWO = 2.0E+0, ZERO = 0.0E+0 )
  151. * ..
  152. * .. Local Scalars ..
  153. LOGICAL UPPER
  154. CHARACTER DIST, PACKIT, TYPE
  155. CHARACTER*3 PATH
  156. INTEGER I, IY, J, JC, JCNEXT, JCOUNT, JJ, JL, JR, JX,
  157. $ KL, KU, MODE
  158. REAL ANORM, BIGNUM, BNORM, BSCAL, C, CNDNUM, REXP,
  159. $ SFAC, SMLNUM, T, TEXP, TLEFT, TSCAL, ULP, UNFL,
  160. $ X, Y, Z
  161. COMPLEX CTEMP, PLUS1, PLUS2, RA, RB, S, STAR1
  162. * ..
  163. * .. External Functions ..
  164. LOGICAL LSAME
  165. INTEGER ICAMAX
  166. REAL SLAMCH
  167. COMPLEX CLARND
  168. EXTERNAL LSAME, ICAMAX, SLAMCH, CLARND
  169. * ..
  170. * .. External Subroutines ..
  171. EXTERNAL CLARNV, CLATB4, CLATMS, CROT, CROTG, CSSCAL,
  172. $ SLABAD, SLARNV
  173. * ..
  174. * .. Intrinsic Functions ..
  175. INTRINSIC ABS, CMPLX, CONJG, MAX, REAL, SQRT
  176. * ..
  177. * .. Executable Statements ..
  178. *
  179. PATH( 1: 1 ) = 'Complex precision'
  180. PATH( 2: 3 ) = 'TP'
  181. UNFL = SLAMCH( 'Safe minimum' )
  182. ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
  183. SMLNUM = UNFL
  184. BIGNUM = ( ONE-ULP ) / SMLNUM
  185. CALL SLABAD( SMLNUM, BIGNUM )
  186. IF( ( IMAT.GE.7 .AND. IMAT.LE.10 ) .OR. IMAT.EQ.18 ) THEN
  187. DIAG = 'U'
  188. ELSE
  189. DIAG = 'N'
  190. END IF
  191. INFO = 0
  192. *
  193. * Quick return if N.LE.0.
  194. *
  195. IF( N.LE.0 )
  196. $ RETURN
  197. *
  198. * Call CLATB4 to set parameters for CLATMS.
  199. *
  200. UPPER = LSAME( UPLO, 'U' )
  201. IF( UPPER ) THEN
  202. CALL CLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
  203. $ CNDNUM, DIST )
  204. PACKIT = 'C'
  205. ELSE
  206. CALL CLATB4( PATH, -IMAT, N, N, TYPE, KL, KU, ANORM, MODE,
  207. $ CNDNUM, DIST )
  208. PACKIT = 'R'
  209. END IF
  210. *
  211. * IMAT <= 6: Non-unit triangular matrix
  212. *
  213. IF( IMAT.LE.6 ) THEN
  214. CALL CLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE, CNDNUM,
  215. $ ANORM, KL, KU, PACKIT, AP, N, WORK, INFO )
  216. *
  217. * IMAT > 6: Unit triangular matrix
  218. * The diagonal is deliberately set to something other than 1.
  219. *
  220. * IMAT = 7: Matrix is the identity
  221. *
  222. ELSE IF( IMAT.EQ.7 ) THEN
  223. IF( UPPER ) THEN
  224. JC = 1
  225. DO 20 J = 1, N
  226. DO 10 I = 1, J - 1
  227. AP( JC+I-1 ) = ZERO
  228. 10 CONTINUE
  229. AP( JC+J-1 ) = J
  230. JC = JC + J
  231. 20 CONTINUE
  232. ELSE
  233. JC = 1
  234. DO 40 J = 1, N
  235. AP( JC ) = J
  236. DO 30 I = J + 1, N
  237. AP( JC+I-J ) = ZERO
  238. 30 CONTINUE
  239. JC = JC + N - J + 1
  240. 40 CONTINUE
  241. END IF
  242. *
  243. * IMAT > 7: Non-trivial unit triangular matrix
  244. *
  245. * Generate a unit triangular matrix T with condition CNDNUM by
  246. * forming a triangular matrix with known singular values and
  247. * filling in the zero entries with Givens rotations.
  248. *
  249. ELSE IF( IMAT.LE.10 ) THEN
  250. IF( UPPER ) THEN
  251. JC = 0
  252. DO 60 J = 1, N
  253. DO 50 I = 1, J - 1
  254. AP( JC+I ) = ZERO
  255. 50 CONTINUE
  256. AP( JC+J ) = J
  257. JC = JC + J
  258. 60 CONTINUE
  259. ELSE
  260. JC = 1
  261. DO 80 J = 1, N
  262. AP( JC ) = J
  263. DO 70 I = J + 1, N
  264. AP( JC+I-J ) = ZERO
  265. 70 CONTINUE
  266. JC = JC + N - J + 1
  267. 80 CONTINUE
  268. END IF
  269. *
  270. * Since the trace of a unit triangular matrix is 1, the product
  271. * of its singular values must be 1. Let s = sqrt(CNDNUM),
  272. * x = sqrt(s) - 1/sqrt(s), y = sqrt(2/(n-2))*x, and z = x**2.
  273. * The following triangular matrix has singular values s, 1, 1,
  274. * ..., 1, 1/s:
  275. *
  276. * 1 y y y ... y y z
  277. * 1 0 0 ... 0 0 y
  278. * 1 0 ... 0 0 y
  279. * . ... . . .
  280. * . . . .
  281. * 1 0 y
  282. * 1 y
  283. * 1
  284. *
  285. * To fill in the zeros, we first multiply by a matrix with small
  286. * condition number of the form
  287. *
  288. * 1 0 0 0 0 ...
  289. * 1 + * 0 0 ...
  290. * 1 + 0 0 0
  291. * 1 + * 0 0
  292. * 1 + 0 0
  293. * ...
  294. * 1 + 0
  295. * 1 0
  296. * 1
  297. *
  298. * Each element marked with a '*' is formed by taking the product
  299. * of the adjacent elements marked with '+'. The '*'s can be
  300. * chosen freely, and the '+'s are chosen so that the inverse of
  301. * T will have elements of the same magnitude as T. If the *'s in
  302. * both T and inv(T) have small magnitude, T is well conditioned.
  303. * The two offdiagonals of T are stored in WORK.
  304. *
  305. * The product of these two matrices has the form
  306. *
  307. * 1 y y y y y . y y z
  308. * 1 + * 0 0 . 0 0 y
  309. * 1 + 0 0 . 0 0 y
  310. * 1 + * . . . .
  311. * 1 + . . . .
  312. * . . . . .
  313. * . . . .
  314. * 1 + y
  315. * 1 y
  316. * 1
  317. *
  318. * Now we multiply by Givens rotations, using the fact that
  319. *
  320. * [ c s ] [ 1 w ] [ -c -s ] = [ 1 -w ]
  321. * [ -s c ] [ 0 1 ] [ s -c ] [ 0 1 ]
  322. * and
  323. * [ -c -s ] [ 1 0 ] [ c s ] = [ 1 0 ]
  324. * [ s -c ] [ w 1 ] [ -s c ] [ -w 1 ]
  325. *
  326. * where c = w / sqrt(w**2+4) and s = 2 / sqrt(w**2+4).
  327. *
  328. STAR1 = 0.25*CLARND( 5, ISEED )
  329. SFAC = 0.5
  330. PLUS1 = SFAC*CLARND( 5, ISEED )
  331. DO 90 J = 1, N, 2
  332. PLUS2 = STAR1 / PLUS1
  333. WORK( J ) = PLUS1
  334. WORK( N+J ) = STAR1
  335. IF( J+1.LE.N ) THEN
  336. WORK( J+1 ) = PLUS2
  337. WORK( N+J+1 ) = ZERO
  338. PLUS1 = STAR1 / PLUS2
  339. REXP = REAL( CLARND( 2, ISEED ) )
  340. IF( REXP.LT.ZERO ) THEN
  341. STAR1 = -SFAC**( ONE-REXP )*CLARND( 5, ISEED )
  342. ELSE
  343. STAR1 = SFAC**( ONE+REXP )*CLARND( 5, ISEED )
  344. END IF
  345. END IF
  346. 90 CONTINUE
  347. *
  348. X = SQRT( CNDNUM ) - ONE / SQRT( CNDNUM )
  349. IF( N.GT.2 ) THEN
  350. Y = SQRT( TWO / REAL( N-2 ) )*X
  351. ELSE
  352. Y = ZERO
  353. END IF
  354. Z = X*X
  355. *
  356. IF( UPPER ) THEN
  357. *
  358. * Set the upper triangle of A with a unit triangular matrix
  359. * of known condition number.
  360. *
  361. JC = 1
  362. DO 100 J = 2, N
  363. AP( JC+1 ) = Y
  364. IF( J.GT.2 )
  365. $ AP( JC+J-1 ) = WORK( J-2 )
  366. IF( J.GT.3 )
  367. $ AP( JC+J-2 ) = WORK( N+J-3 )
  368. JC = JC + J
  369. 100 CONTINUE
  370. JC = JC - N
  371. AP( JC+1 ) = Z
  372. DO 110 J = 2, N - 1
  373. AP( JC+J ) = Y
  374. 110 CONTINUE
  375. ELSE
  376. *
  377. * Set the lower triangle of A with a unit triangular matrix
  378. * of known condition number.
  379. *
  380. DO 120 I = 2, N - 1
  381. AP( I ) = Y
  382. 120 CONTINUE
  383. AP( N ) = Z
  384. JC = N + 1
  385. DO 130 J = 2, N - 1
  386. AP( JC+1 ) = WORK( J-1 )
  387. IF( J.LT.N-1 )
  388. $ AP( JC+2 ) = WORK( N+J-1 )
  389. AP( JC+N-J ) = Y
  390. JC = JC + N - J + 1
  391. 130 CONTINUE
  392. END IF
  393. *
  394. * Fill in the zeros using Givens rotations
  395. *
  396. IF( UPPER ) THEN
  397. JC = 1
  398. DO 150 J = 1, N - 1
  399. JCNEXT = JC + J
  400. RA = AP( JCNEXT+J-1 )
  401. RB = TWO
  402. CALL CROTG( RA, RB, C, S )
  403. *
  404. * Multiply by [ c s; -conjg(s) c] on the left.
  405. *
  406. IF( N.GT.J+1 ) THEN
  407. JX = JCNEXT + J
  408. DO 140 I = J + 2, N
  409. CTEMP = C*AP( JX+J ) + S*AP( JX+J+1 )
  410. AP( JX+J+1 ) = -CONJG( S )*AP( JX+J ) +
  411. $ C*AP( JX+J+1 )
  412. AP( JX+J ) = CTEMP
  413. JX = JX + I
  414. 140 CONTINUE
  415. END IF
  416. *
  417. * Multiply by [-c -s; conjg(s) -c] on the right.
  418. *
  419. IF( J.GT.1 )
  420. $ CALL CROT( J-1, AP( JCNEXT ), 1, AP( JC ), 1, -C, -S )
  421. *
  422. * Negate A(J,J+1).
  423. *
  424. AP( JCNEXT+J-1 ) = -AP( JCNEXT+J-1 )
  425. JC = JCNEXT
  426. 150 CONTINUE
  427. ELSE
  428. JC = 1
  429. DO 170 J = 1, N - 1
  430. JCNEXT = JC + N - J + 1
  431. RA = AP( JC+1 )
  432. RB = TWO
  433. CALL CROTG( RA, RB, C, S )
  434. S = CONJG( S )
  435. *
  436. * Multiply by [ c -s; conjg(s) c] on the right.
  437. *
  438. IF( N.GT.J+1 )
  439. $ CALL CROT( N-J-1, AP( JCNEXT+1 ), 1, AP( JC+2 ), 1, C,
  440. $ -S )
  441. *
  442. * Multiply by [-c s; -conjg(s) -c] on the left.
  443. *
  444. IF( J.GT.1 ) THEN
  445. JX = 1
  446. DO 160 I = 1, J - 1
  447. CTEMP = -C*AP( JX+J-I ) + S*AP( JX+J-I+1 )
  448. AP( JX+J-I+1 ) = -CONJG( S )*AP( JX+J-I ) -
  449. $ C*AP( JX+J-I+1 )
  450. AP( JX+J-I ) = CTEMP
  451. JX = JX + N - I + 1
  452. 160 CONTINUE
  453. END IF
  454. *
  455. * Negate A(J+1,J).
  456. *
  457. AP( JC+1 ) = -AP( JC+1 )
  458. JC = JCNEXT
  459. 170 CONTINUE
  460. END IF
  461. *
  462. * IMAT > 10: Pathological test cases. These triangular matrices
  463. * are badly scaled or badly conditioned, so when used in solving a
  464. * triangular system they may cause overflow in the solution vector.
  465. *
  466. ELSE IF( IMAT.EQ.11 ) THEN
  467. *
  468. * Type 11: Generate a triangular matrix with elements between
  469. * -1 and 1. Give the diagonal norm 2 to make it well-conditioned.
  470. * Make the right hand side large so that it requires scaling.
  471. *
  472. IF( UPPER ) THEN
  473. JC = 1
  474. DO 180 J = 1, N
  475. CALL CLARNV( 4, ISEED, J-1, AP( JC ) )
  476. AP( JC+J-1 ) = CLARND( 5, ISEED )*TWO
  477. JC = JC + J
  478. 180 CONTINUE
  479. ELSE
  480. JC = 1
  481. DO 190 J = 1, N
  482. IF( J.LT.N )
  483. $ CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) )
  484. AP( JC ) = CLARND( 5, ISEED )*TWO
  485. JC = JC + N - J + 1
  486. 190 CONTINUE
  487. END IF
  488. *
  489. * Set the right hand side so that the largest value is BIGNUM.
  490. *
  491. CALL CLARNV( 2, ISEED, N, B )
  492. IY = ICAMAX( N, B, 1 )
  493. BNORM = ABS( B( IY ) )
  494. BSCAL = BIGNUM / MAX( ONE, BNORM )
  495. CALL CSSCAL( N, BSCAL, B, 1 )
  496. *
  497. ELSE IF( IMAT.EQ.12 ) THEN
  498. *
  499. * Type 12: Make the first diagonal element in the solve small to
  500. * cause immediate overflow when dividing by T(j,j).
  501. * In type 12, the offdiagonal elements are small (CNORM(j) < 1).
  502. *
  503. CALL CLARNV( 2, ISEED, N, B )
  504. TSCAL = ONE / MAX( ONE, REAL( N-1 ) )
  505. IF( UPPER ) THEN
  506. JC = 1
  507. DO 200 J = 1, N
  508. CALL CLARNV( 4, ISEED, J-1, AP( JC ) )
  509. CALL CSSCAL( J-1, TSCAL, AP( JC ), 1 )
  510. AP( JC+J-1 ) = CLARND( 5, ISEED )
  511. JC = JC + J
  512. 200 CONTINUE
  513. AP( N*( N+1 ) / 2 ) = SMLNUM*AP( N*( N+1 ) / 2 )
  514. ELSE
  515. JC = 1
  516. DO 210 J = 1, N
  517. CALL CLARNV( 2, ISEED, N-J, AP( JC+1 ) )
  518. CALL CSSCAL( N-J, TSCAL, AP( JC+1 ), 1 )
  519. AP( JC ) = CLARND( 5, ISEED )
  520. JC = JC + N - J + 1
  521. 210 CONTINUE
  522. AP( 1 ) = SMLNUM*AP( 1 )
  523. END IF
  524. *
  525. ELSE IF( IMAT.EQ.13 ) THEN
  526. *
  527. * Type 13: Make the first diagonal element in the solve small to
  528. * cause immediate overflow when dividing by T(j,j).
  529. * In type 13, the offdiagonal elements are O(1) (CNORM(j) > 1).
  530. *
  531. CALL CLARNV( 2, ISEED, N, B )
  532. IF( UPPER ) THEN
  533. JC = 1
  534. DO 220 J = 1, N
  535. CALL CLARNV( 4, ISEED, J-1, AP( JC ) )
  536. AP( JC+J-1 ) = CLARND( 5, ISEED )
  537. JC = JC + J
  538. 220 CONTINUE
  539. AP( N*( N+1 ) / 2 ) = SMLNUM*AP( N*( N+1 ) / 2 )
  540. ELSE
  541. JC = 1
  542. DO 230 J = 1, N
  543. CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) )
  544. AP( JC ) = CLARND( 5, ISEED )
  545. JC = JC + N - J + 1
  546. 230 CONTINUE
  547. AP( 1 ) = SMLNUM*AP( 1 )
  548. END IF
  549. *
  550. ELSE IF( IMAT.EQ.14 ) THEN
  551. *
  552. * Type 14: T is diagonal with small numbers on the diagonal to
  553. * make the growth factor underflow, but a small right hand side
  554. * chosen so that the solution does not overflow.
  555. *
  556. IF( UPPER ) THEN
  557. JCOUNT = 1
  558. JC = ( N-1 )*N / 2 + 1
  559. DO 250 J = N, 1, -1
  560. DO 240 I = 1, J - 1
  561. AP( JC+I-1 ) = ZERO
  562. 240 CONTINUE
  563. IF( JCOUNT.LE.2 ) THEN
  564. AP( JC+J-1 ) = SMLNUM*CLARND( 5, ISEED )
  565. ELSE
  566. AP( JC+J-1 ) = CLARND( 5, ISEED )
  567. END IF
  568. JCOUNT = JCOUNT + 1
  569. IF( JCOUNT.GT.4 )
  570. $ JCOUNT = 1
  571. JC = JC - J + 1
  572. 250 CONTINUE
  573. ELSE
  574. JCOUNT = 1
  575. JC = 1
  576. DO 270 J = 1, N
  577. DO 260 I = J + 1, N
  578. AP( JC+I-J ) = ZERO
  579. 260 CONTINUE
  580. IF( JCOUNT.LE.2 ) THEN
  581. AP( JC ) = SMLNUM*CLARND( 5, ISEED )
  582. ELSE
  583. AP( JC ) = CLARND( 5, ISEED )
  584. END IF
  585. JCOUNT = JCOUNT + 1
  586. IF( JCOUNT.GT.4 )
  587. $ JCOUNT = 1
  588. JC = JC + N - J + 1
  589. 270 CONTINUE
  590. END IF
  591. *
  592. * Set the right hand side alternately zero and small.
  593. *
  594. IF( UPPER ) THEN
  595. B( 1 ) = ZERO
  596. DO 280 I = N, 2, -2
  597. B( I ) = ZERO
  598. B( I-1 ) = SMLNUM*CLARND( 5, ISEED )
  599. 280 CONTINUE
  600. ELSE
  601. B( N ) = ZERO
  602. DO 290 I = 1, N - 1, 2
  603. B( I ) = ZERO
  604. B( I+1 ) = SMLNUM*CLARND( 5, ISEED )
  605. 290 CONTINUE
  606. END IF
  607. *
  608. ELSE IF( IMAT.EQ.15 ) THEN
  609. *
  610. * Type 15: Make the diagonal elements small to cause gradual
  611. * overflow when dividing by T(j,j). To control the amount of
  612. * scaling needed, the matrix is bidiagonal.
  613. *
  614. TEXP = ONE / MAX( ONE, REAL( N-1 ) )
  615. TSCAL = SMLNUM**TEXP
  616. CALL CLARNV( 4, ISEED, N, B )
  617. IF( UPPER ) THEN
  618. JC = 1
  619. DO 310 J = 1, N
  620. DO 300 I = 1, J - 2
  621. AP( JC+I-1 ) = ZERO
  622. 300 CONTINUE
  623. IF( J.GT.1 )
  624. $ AP( JC+J-2 ) = CMPLX( -ONE, -ONE )
  625. AP( JC+J-1 ) = TSCAL*CLARND( 5, ISEED )
  626. JC = JC + J
  627. 310 CONTINUE
  628. B( N ) = CMPLX( ONE, ONE )
  629. ELSE
  630. JC = 1
  631. DO 330 J = 1, N
  632. DO 320 I = J + 2, N
  633. AP( JC+I-J ) = ZERO
  634. 320 CONTINUE
  635. IF( J.LT.N )
  636. $ AP( JC+1 ) = CMPLX( -ONE, -ONE )
  637. AP( JC ) = TSCAL*CLARND( 5, ISEED )
  638. JC = JC + N - J + 1
  639. 330 CONTINUE
  640. B( 1 ) = CMPLX( ONE, ONE )
  641. END IF
  642. *
  643. ELSE IF( IMAT.EQ.16 ) THEN
  644. *
  645. * Type 16: One zero diagonal element.
  646. *
  647. IY = N / 2 + 1
  648. IF( UPPER ) THEN
  649. JC = 1
  650. DO 340 J = 1, N
  651. CALL CLARNV( 4, ISEED, J, AP( JC ) )
  652. IF( J.NE.IY ) THEN
  653. AP( JC+J-1 ) = CLARND( 5, ISEED )*TWO
  654. ELSE
  655. AP( JC+J-1 ) = ZERO
  656. END IF
  657. JC = JC + J
  658. 340 CONTINUE
  659. ELSE
  660. JC = 1
  661. DO 350 J = 1, N
  662. CALL CLARNV( 4, ISEED, N-J+1, AP( JC ) )
  663. IF( J.NE.IY ) THEN
  664. AP( JC ) = CLARND( 5, ISEED )*TWO
  665. ELSE
  666. AP( JC ) = ZERO
  667. END IF
  668. JC = JC + N - J + 1
  669. 350 CONTINUE
  670. END IF
  671. CALL CLARNV( 2, ISEED, N, B )
  672. CALL CSSCAL( N, TWO, B, 1 )
  673. *
  674. ELSE IF( IMAT.EQ.17 ) THEN
  675. *
  676. * Type 17: Make the offdiagonal elements large to cause overflow
  677. * when adding a column of T. In the non-transposed case, the
  678. * matrix is constructed to cause overflow when adding a column in
  679. * every other step.
  680. *
  681. TSCAL = UNFL / ULP
  682. TSCAL = ( ONE-ULP ) / TSCAL
  683. DO 360 J = 1, N*( N+1 ) / 2
  684. AP( J ) = ZERO
  685. 360 CONTINUE
  686. TEXP = ONE
  687. IF( UPPER ) THEN
  688. JC = ( N-1 )*N / 2 + 1
  689. DO 370 J = N, 2, -2
  690. AP( JC ) = -TSCAL / REAL( N+1 )
  691. AP( JC+J-1 ) = ONE
  692. B( J ) = TEXP*( ONE-ULP )
  693. JC = JC - J + 1
  694. AP( JC ) = -( TSCAL / REAL( N+1 ) ) / REAL( N+2 )
  695. AP( JC+J-2 ) = ONE
  696. B( J-1 ) = TEXP*REAL( N*N+N-1 )
  697. TEXP = TEXP*TWO
  698. JC = JC - J + 2
  699. 370 CONTINUE
  700. B( 1 ) = ( REAL( N+1 ) / REAL( N+2 ) )*TSCAL
  701. ELSE
  702. JC = 1
  703. DO 380 J = 1, N - 1, 2
  704. AP( JC+N-J ) = -TSCAL / REAL( N+1 )
  705. AP( JC ) = ONE
  706. B( J ) = TEXP*( ONE-ULP )
  707. JC = JC + N - J + 1
  708. AP( JC+N-J-1 ) = -( TSCAL / REAL( N+1 ) ) / REAL( N+2 )
  709. AP( JC ) = ONE
  710. B( J+1 ) = TEXP*REAL( N*N+N-1 )
  711. TEXP = TEXP*TWO
  712. JC = JC + N - J
  713. 380 CONTINUE
  714. B( N ) = ( REAL( N+1 ) / REAL( N+2 ) )*TSCAL
  715. END IF
  716. *
  717. ELSE IF( IMAT.EQ.18 ) THEN
  718. *
  719. * Type 18: Generate a unit triangular matrix with elements
  720. * between -1 and 1, and make the right hand side large so that it
  721. * requires scaling.
  722. *
  723. IF( UPPER ) THEN
  724. JC = 1
  725. DO 390 J = 1, N
  726. CALL CLARNV( 4, ISEED, J-1, AP( JC ) )
  727. AP( JC+J-1 ) = ZERO
  728. JC = JC + J
  729. 390 CONTINUE
  730. ELSE
  731. JC = 1
  732. DO 400 J = 1, N
  733. IF( J.LT.N )
  734. $ CALL CLARNV( 4, ISEED, N-J, AP( JC+1 ) )
  735. AP( JC ) = ZERO
  736. JC = JC + N - J + 1
  737. 400 CONTINUE
  738. END IF
  739. *
  740. * Set the right hand side so that the largest value is BIGNUM.
  741. *
  742. CALL CLARNV( 2, ISEED, N, B )
  743. IY = ICAMAX( N, B, 1 )
  744. BNORM = ABS( B( IY ) )
  745. BSCAL = BIGNUM / MAX( ONE, BNORM )
  746. CALL CSSCAL( N, BSCAL, B, 1 )
  747. *
  748. ELSE IF( IMAT.EQ.19 ) THEN
  749. *
  750. * Type 19: Generate a triangular matrix with elements between
  751. * BIGNUM/(n-1) and BIGNUM so that at least one of the column
  752. * norms will exceed BIGNUM.
  753. * 1/3/91: CLATPS no longer can handle this case
  754. *
  755. TLEFT = BIGNUM / MAX( ONE, REAL( N-1 ) )
  756. TSCAL = BIGNUM*( REAL( N-1 ) / MAX( ONE, REAL( N ) ) )
  757. IF( UPPER ) THEN
  758. JC = 1
  759. DO 420 J = 1, N
  760. CALL CLARNV( 5, ISEED, J, AP( JC ) )
  761. CALL SLARNV( 1, ISEED, J, RWORK )
  762. DO 410 I = 1, J
  763. AP( JC+I-1 ) = AP( JC+I-1 )*( TLEFT+RWORK( I )*TSCAL )
  764. 410 CONTINUE
  765. JC = JC + J
  766. 420 CONTINUE
  767. ELSE
  768. JC = 1
  769. DO 440 J = 1, N
  770. CALL CLARNV( 5, ISEED, N-J+1, AP( JC ) )
  771. CALL SLARNV( 1, ISEED, N-J+1, RWORK )
  772. DO 430 I = J, N
  773. AP( JC+I-J ) = AP( JC+I-J )*
  774. $ ( TLEFT+RWORK( I-J+1 )*TSCAL )
  775. 430 CONTINUE
  776. JC = JC + N - J + 1
  777. 440 CONTINUE
  778. END IF
  779. CALL CLARNV( 2, ISEED, N, B )
  780. CALL CSSCAL( N, TWO, B, 1 )
  781. END IF
  782. *
  783. * Flip the matrix across its counter-diagonal if the transpose will
  784. * be used.
  785. *
  786. IF( .NOT.LSAME( TRANS, 'N' ) ) THEN
  787. IF( UPPER ) THEN
  788. JJ = 1
  789. JR = N*( N+1 ) / 2
  790. DO 460 J = 1, N / 2
  791. JL = JJ
  792. DO 450 I = J, N - J
  793. T = REAL( AP( JR-I+J ) )
  794. AP( JR-I+J ) = AP( JL )
  795. AP( JL ) = T
  796. JL = JL + I
  797. 450 CONTINUE
  798. JJ = JJ + J + 1
  799. JR = JR - ( N-J+1 )
  800. 460 CONTINUE
  801. ELSE
  802. JL = 1
  803. JJ = N*( N+1 ) / 2
  804. DO 480 J = 1, N / 2
  805. JR = JJ
  806. DO 470 I = J, N - J
  807. T = REAL( AP( JL+I-J ) )
  808. AP( JL+I-J ) = AP( JR )
  809. AP( JR ) = T
  810. JR = JR - I
  811. 470 CONTINUE
  812. JL = JL + N - J + 1
  813. JJ = JJ - J - 1
  814. 480 CONTINUE
  815. END IF
  816. END IF
  817. *
  818. RETURN
  819. *
  820. * End of CLATTP
  821. *
  822. END