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