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.

dlattp.f 25 kB

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