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