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