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