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