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 26 kB

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