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