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