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