LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlatmt.f
Go to the documentation of this file.
1*> \brief \b DLATMT
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 DLATMT( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
12* RANK, KL, KU, PACK, A, LDA, WORK, INFO )
13*
14* .. Scalar Arguments ..
15* DOUBLE PRECISION COND, DMAX
16* INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK
17* CHARACTER DIST, PACK, SYM
18* ..
19* .. Array Arguments ..
20* DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
21* INTEGER ISEED( 4 )
22* ..
23*
24*
25*> \par Purpose:
26* =============
27*>
28*> \verbatim
29*>
30*> DLATMT generates random matrices with specified singular values
31*> (or symmetric/hermitian with specified eigenvalues)
32*> for testing LAPACK programs.
33*>
34*> DLATMT operates by applying the following sequence of
35*> operations:
36*>
37*> Set the diagonal to D, where D may be input or
38*> computed according to MODE, COND, DMAX, and SYM
39*> as described below.
40*>
41*> Generate a matrix with the appropriate band structure, by one
42*> of two methods:
43*>
44*> Method A:
45*> Generate a dense M x N matrix by multiplying D on the left
46*> and the right by random unitary matrices, then:
47*>
48*> Reduce the bandwidth according to KL and KU, using
49*> Householder transformations.
50*>
51*> Method B:
52*> Convert the bandwidth-0 (i.e., diagonal) matrix to a
53*> bandwidth-1 matrix using Givens rotations, "chasing"
54*> out-of-band elements back, much as in QR; then
55*> convert the bandwidth-1 to a bandwidth-2 matrix, etc.
56*> Note that for reasonably small bandwidths (relative to
57*> M and N) this requires less storage, as a dense matrix
58*> is not generated. Also, for symmetric matrices, only
59*> one triangle is generated.
60*>
61*> Method A is chosen if the bandwidth is a large fraction of the
62*> order of the matrix, and LDA is at least M (so a dense
63*> matrix can be stored.) Method B is chosen if the bandwidth
64*> is small (< 1/2 N for symmetric, < .3 N+M for
65*> non-symmetric), or LDA is less than M and not less than the
66*> bandwidth.
67*>
68*> Pack the matrix if desired. Options specified by PACK are:
69*> no packing
70*> zero out upper half (if symmetric)
71*> zero out lower half (if symmetric)
72*> store the upper half columnwise (if symmetric or upper
73*> triangular)
74*> store the lower half columnwise (if symmetric or lower
75*> triangular)
76*> store the lower triangle in banded format (if symmetric
77*> or lower triangular)
78*> store the upper triangle in banded format (if symmetric
79*> or upper triangular)
80*> store the entire matrix in banded format
81*> If Method B is chosen, and band format is specified, then the
82*> matrix will be generated in the band format, so no repacking
83*> will be necessary.
84*> \endverbatim
85*
86* Arguments:
87* ==========
88*
89*> \param[in] M
90*> \verbatim
91*> M is INTEGER
92*> The number of rows of A. Not modified.
93*> \endverbatim
94*>
95*> \param[in] N
96*> \verbatim
97*> N is INTEGER
98*> The number of columns of A. Not modified.
99*> \endverbatim
100*>
101*> \param[in] DIST
102*> \verbatim
103*> DIST is CHARACTER*1
104*> On entry, DIST specifies the type of distribution to be used
105*> to generate the random eigen-/singular values.
106*> 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
107*> 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
108*> 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
109*> Not modified.
110*> \endverbatim
111*>
112*> \param[in,out] ISEED
113*> \verbatim
114*> ISEED is INTEGER array, dimension ( 4 )
115*> On entry ISEED specifies the seed of the random number
116*> generator. They should lie between 0 and 4095 inclusive,
117*> and ISEED(4) should be odd. The random number generator
118*> uses a linear congruential sequence limited to small
119*> integers, and so should produce machine independent
120*> random numbers. The values of ISEED are changed on
121*> exit, and can be used in the next call to DLATMT
122*> to continue the same random number sequence.
123*> Changed on exit.
124*> \endverbatim
125*>
126*> \param[in] SYM
127*> \verbatim
128*> SYM is CHARACTER*1
129*> If SYM='S' or 'H', the generated matrix is symmetric, with
130*> eigenvalues specified by D, COND, MODE, and DMAX; they
131*> may be positive, negative, or zero.
132*> If SYM='P', the generated matrix is symmetric, with
133*> eigenvalues (= singular values) specified by D, COND,
134*> MODE, and DMAX; they will not be negative.
135*> If SYM='N', the generated matrix is nonsymmetric, with
136*> singular values specified by D, COND, MODE, and DMAX;
137*> they will not be negative.
138*> Not modified.
139*> \endverbatim
140*>
141*> \param[in,out] D
142*> \verbatim
143*> D is DOUBLE PRECISION array, dimension ( MIN( M , N ) )
144*> This array is used to specify the singular values or
145*> eigenvalues of A (see SYM, above.) If MODE=0, then D is
146*> assumed to contain the singular/eigenvalues, otherwise
147*> they will be computed according to MODE, COND, and DMAX,
148*> and placed in D.
149*> Modified if MODE is nonzero.
150*> \endverbatim
151*>
152*> \param[in] MODE
153*> \verbatim
154*> MODE is INTEGER
155*> On entry this describes how the singular/eigenvalues are to
156*> be specified:
157*> MODE = 0 means use D as input
158*>
159*> MODE = 1 sets D(1)=1 and D(2:RANK)=1.0/COND
160*> MODE = 2 sets D(1:RANK-1)=1 and D(RANK)=1.0/COND
161*> MODE = 3 sets D(I)=COND**(-(I-1)/(RANK-1))
162*>
163*> MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
164*> MODE = 5 sets D to random numbers in the range
165*> ( 1/COND , 1 ) such that their logarithms
166*> are uniformly distributed.
167*> MODE = 6 set D to random numbers from same distribution
168*> as the rest of the matrix.
169*> MODE < 0 has the same meaning as ABS(MODE), except that
170*> the order of the elements of D is reversed.
171*> Thus if MODE is positive, D has entries ranging from
172*> 1 to 1/COND, if negative, from 1/COND to 1,
173*> If SYM='S' or 'H', and MODE is neither 0, 6, nor -6, then
174*> the elements of D will also be multiplied by a random
175*> sign (i.e., +1 or -1.)
176*> Not modified.
177*> \endverbatim
178*>
179*> \param[in] COND
180*> \verbatim
181*> COND is DOUBLE PRECISION
182*> On entry, this is used as described under MODE above.
183*> If used, it must be >= 1. Not modified.
184*> \endverbatim
185*>
186*> \param[in] DMAX
187*> \verbatim
188*> DMAX is DOUBLE PRECISION
189*> If MODE is neither -6, 0 nor 6, the contents of D, as
190*> computed according to MODE and COND, will be scaled by
191*> DMAX / max(abs(D(i))); thus, the maximum absolute eigen- or
192*> singular value (which is to say the norm) will be abs(DMAX).
193*> Note that DMAX need not be positive: if DMAX is negative
194*> (or zero), D will be scaled by a negative number (or zero).
195*> Not modified.
196*> \endverbatim
197*>
198*> \param[in] RANK
199*> \verbatim
200*> RANK is INTEGER
201*> The rank of matrix to be generated for modes 1,2,3 only.
202*> D( RANK+1:N ) = 0.
203*> Not modified.
204*> \endverbatim
205*>
206*> \param[in] KL
207*> \verbatim
208*> KL is INTEGER
209*> This specifies the lower bandwidth of the matrix. For
210*> example, KL=0 implies upper triangular, KL=1 implies upper
211*> Hessenberg, and KL being at least M-1 means that the matrix
212*> has full lower bandwidth. KL must equal KU if the matrix
213*> is symmetric.
214*> Not modified.
215*> \endverbatim
216*>
217*> \param[in] KU
218*> \verbatim
219*> KU is INTEGER
220*> This specifies the upper bandwidth of the matrix. For
221*> example, KU=0 implies lower triangular, KU=1 implies lower
222*> Hessenberg, and KU being at least N-1 means that the matrix
223*> has full upper bandwidth. KL must equal KU if the matrix
224*> is symmetric.
225*> Not modified.
226*> \endverbatim
227*>
228*> \param[in] PACK
229*> \verbatim
230*> PACK is CHARACTER*1
231*> This specifies packing of matrix as follows:
232*> 'N' => no packing
233*> 'U' => zero out all subdiagonal entries (if symmetric)
234*> 'L' => zero out all superdiagonal entries (if symmetric)
235*> 'C' => store the upper triangle columnwise
236*> (only if the matrix is symmetric or upper triangular)
237*> 'R' => store the lower triangle columnwise
238*> (only if the matrix is symmetric or lower triangular)
239*> 'B' => store the lower triangle in band storage scheme
240*> (only if matrix symmetric or lower triangular)
241*> 'Q' => store the upper triangle in band storage scheme
242*> (only if matrix symmetric or upper triangular)
243*> 'Z' => store the entire matrix in band storage scheme
244*> (pivoting can be provided for by using this
245*> option to store A in the trailing rows of
246*> the allocated storage)
247*>
248*> Using these options, the various LAPACK packed and banded
249*> storage schemes can be obtained:
250*> GB - use 'Z'
251*> PB, SB or TB - use 'B' or 'Q'
252*> PP, SP or TP - use 'C' or 'R'
253*>
254*> If two calls to DLATMT differ only in the PACK parameter,
255*> they will generate mathematically equivalent matrices.
256*> Not modified.
257*> \endverbatim
258*>
259*> \param[in,out] A
260*> \verbatim
261*> A is DOUBLE PRECISION array, dimension ( LDA, N )
262*> On exit A is the desired test matrix. A is first generated
263*> in full (unpacked) form, and then packed, if so specified
264*> by PACK. Thus, the first M elements of the first N
265*> columns will always be modified. If PACK specifies a
266*> packed or banded storage scheme, all LDA elements of the
267*> first N columns will be modified; the elements of the
268*> array which do not correspond to elements of the generated
269*> matrix are set to zero.
270*> Modified.
271*> \endverbatim
272*>
273*> \param[in] LDA
274*> \verbatim
275*> LDA is INTEGER
276*> LDA specifies the first dimension of A as declared in the
277*> calling program. If PACK='N', 'U', 'L', 'C', or 'R', then
278*> LDA must be at least M. If PACK='B' or 'Q', then LDA must
279*> be at least MIN( KL, M-1) (which is equal to MIN(KU,N-1)).
280*> If PACK='Z', LDA must be large enough to hold the packed
281*> array: MIN( KU, N-1) + MIN( KL, M-1) + 1.
282*> Not modified.
283*> \endverbatim
284*>
285*> \param[out] WORK
286*> \verbatim
287*> WORK is DOUBLE PRECISION array, dimension ( 3*MAX( N , M ) )
288*> Workspace.
289*> Modified.
290*> \endverbatim
291*>
292*> \param[out] INFO
293*> \verbatim
294*> INFO is INTEGER
295*> Error code. On exit, INFO will be set to one of the
296*> following values:
297*> 0 => normal return
298*> -1 => M negative or unequal to N and SYM='S', 'H', or 'P'
299*> -2 => N negative
300*> -3 => DIST illegal string
301*> -5 => SYM illegal string
302*> -7 => MODE not in range -6 to 6
303*> -8 => COND less than 1.0, and MODE neither -6, 0 nor 6
304*> -10 => KL negative
305*> -11 => KU negative, or SYM='S' or 'H' and KU not equal to KL
306*> -12 => PACK illegal string, or PACK='U' or 'L', and SYM='N';
307*> or PACK='C' or 'Q' and SYM='N' and KL is not zero;
308*> or PACK='R' or 'B' and SYM='N' and KU is not zero;
309*> or PACK='U', 'L', 'C', 'R', 'B', or 'Q', and M is not
310*> N.
311*> -14 => LDA is less than M, or PACK='Z' and LDA is less than
312*> MIN(KU,N-1) + MIN(KL,M-1) + 1.
313*> 1 => Error return from DLATM7
314*> 2 => Cannot scale to DMAX (max. sing. value is 0)
315*> 3 => Error return from DLAGGE or DLAGSY
316*> \endverbatim
317*
318* Authors:
319* ========
320*
321*> \author Univ. of Tennessee
322*> \author Univ. of California Berkeley
323*> \author Univ. of Colorado Denver
324*> \author NAG Ltd.
325*
326*> \ingroup double_matgen
327*
328* =====================================================================
329 SUBROUTINE dlatmt( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
330 $ RANK, KL, KU, PACK, A, LDA, WORK, INFO )
331*
332* -- LAPACK computational routine --
333* -- LAPACK is a software package provided by Univ. of Tennessee, --
334* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
335*
336* .. Scalar Arguments ..
337 DOUBLE PRECISION COND, DMAX
338 INTEGER INFO, KL, KU, LDA, M, MODE, N, RANK
339 CHARACTER DIST, PACK, SYM
340* ..
341* .. Array Arguments ..
342 DOUBLE PRECISION A( LDA, * ), D( * ), WORK( * )
343 INTEGER ISEED( 4 )
344* ..
345*
346* =====================================================================
347*
348* .. Parameters ..
349 DOUBLE PRECISION ZERO
350 parameter( zero = 0.0d0 )
351 DOUBLE PRECISION ONE
352 parameter( one = 1.0d0 )
353 DOUBLE PRECISION TWOPI
354 parameter( twopi = 6.28318530717958647692528676655900576839d+0 )
355* ..
356* .. Local Scalars ..
357 DOUBLE PRECISION ALPHA, ANGLE, C, DUMMY, EXTRA, S, TEMP
358 INTEGER I, IC, ICOL, IDIST, IENDCH, IINFO, IL, ILDA,
359 $ ioffg, ioffst, ipack, ipackg, ir, ir1, ir2,
360 $ irow, irsign, iskew, isym, isympk, j, jc, jch,
361 $ jkl, jku, jr, k, llb, minlda, mnmin, mr, nc,
362 $ uub
363 LOGICAL GIVENS, ILEXTR, ILTEMP, TOPDWN
364* ..
365* .. External Functions ..
366 DOUBLE PRECISION DLARND
367 LOGICAL LSAME
368 EXTERNAL dlarnd, lsame
369* ..
370* .. External Subroutines ..
371 EXTERNAL dlatm7, dcopy, dlagge, dlagsy, dlarot,
373* ..
374* .. Intrinsic Functions ..
375 INTRINSIC abs, cos, dble, max, min, mod, sin
376* ..
377* .. Executable Statements ..
378*
379* 1) Decode and Test the input parameters.
380* Initialize flags & seed.
381*
382 info = 0
383*
384* Quick return if possible
385*
386 IF( m.EQ.0 .OR. n.EQ.0 )
387 $ RETURN
388*
389* Decode DIST
390*
391 IF( lsame( dist, 'U' ) ) THEN
392 idist = 1
393 ELSE IF( lsame( dist, 'S' ) ) THEN
394 idist = 2
395 ELSE IF( lsame( dist, 'N' ) ) THEN
396 idist = 3
397 ELSE
398 idist = -1
399 END IF
400*
401* Decode SYM
402*
403 IF( lsame( sym, 'N' ) ) THEN
404 isym = 1
405 irsign = 0
406 ELSE IF( lsame( sym, 'P' ) ) THEN
407 isym = 2
408 irsign = 0
409 ELSE IF( lsame( sym, 'S' ) ) THEN
410 isym = 2
411 irsign = 1
412 ELSE IF( lsame( sym, 'H' ) ) THEN
413 isym = 2
414 irsign = 1
415 ELSE
416 isym = -1
417 END IF
418*
419* Decode PACK
420*
421 isympk = 0
422 IF( lsame( pack, 'N' ) ) THEN
423 ipack = 0
424 ELSE IF( lsame( pack, 'U' ) ) THEN
425 ipack = 1
426 isympk = 1
427 ELSE IF( lsame( pack, 'L' ) ) THEN
428 ipack = 2
429 isympk = 1
430 ELSE IF( lsame( pack, 'C' ) ) THEN
431 ipack = 3
432 isympk = 2
433 ELSE IF( lsame( pack, 'R' ) ) THEN
434 ipack = 4
435 isympk = 3
436 ELSE IF( lsame( pack, 'B' ) ) THEN
437 ipack = 5
438 isympk = 3
439 ELSE IF( lsame( pack, 'Q' ) ) THEN
440 ipack = 6
441 isympk = 2
442 ELSE IF( lsame( pack, 'Z' ) ) THEN
443 ipack = 7
444 ELSE
445 ipack = -1
446 END IF
447*
448* Set certain internal parameters
449*
450 mnmin = min( m, n )
451 llb = min( kl, m-1 )
452 uub = min( ku, n-1 )
453 mr = min( m, n+llb )
454 nc = min( n, m+uub )
455*
456 IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
457 minlda = uub + 1
458 ELSE IF( ipack.EQ.7 ) THEN
459 minlda = llb + uub + 1
460 ELSE
461 minlda = m
462 END IF
463*
464* Use Givens rotation method if bandwidth small enough,
465* or if LDA is too small to store the matrix unpacked.
466*
467 givens = .false.
468 IF( isym.EQ.1 ) THEN
469 IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
470 $ givens = .true.
471 ELSE
472 IF( 2*llb.LT.m )
473 $ givens = .true.
474 END IF
475 IF( lda.LT.m .AND. lda.GE.minlda )
476 $ givens = .true.
477*
478* Set INFO if an error
479*
480 IF( m.LT.0 ) THEN
481 info = -1
482 ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
483 info = -1
484 ELSE IF( n.LT.0 ) THEN
485 info = -2
486 ELSE IF( idist.EQ.-1 ) THEN
487 info = -3
488 ELSE IF( isym.EQ.-1 ) THEN
489 info = -5
490 ELSE IF( abs( mode ).GT.6 ) THEN
491 info = -7
492 ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
493 $ THEN
494 info = -8
495 ELSE IF( kl.LT.0 ) THEN
496 info = -10
497 ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
498 info = -11
499 ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
500 $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
501 $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
502 $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
503 info = -12
504 ELSE IF( lda.LT.max( 1, minlda ) ) THEN
505 info = -14
506 END IF
507*
508 IF( info.NE.0 ) THEN
509 CALL xerbla( 'DLATMT', -info )
510 RETURN
511 END IF
512*
513* Initialize random number generator
514*
515 DO 100 i = 1, 4
516 iseed( i ) = mod( abs( iseed( i ) ), 4096 )
517 100 CONTINUE
518*
519 IF( mod( iseed( 4 ), 2 ).NE.1 )
520 $ iseed( 4 ) = iseed( 4 ) + 1
521*
522* 2) Set up D if indicated.
523*
524* Compute D according to COND and MODE
525*
526 CALL dlatm7( mode, cond, irsign, idist, iseed, d, mnmin, rank,
527 $ iinfo )
528 IF( iinfo.NE.0 ) THEN
529 info = 1
530 RETURN
531 END IF
532*
533* Choose Top-Down if D is (apparently) increasing,
534* Bottom-Up if D is (apparently) decreasing.
535*
536 IF( abs( d( 1 ) ).LE.abs( d( rank ) ) ) THEN
537 topdwn = .true.
538 ELSE
539 topdwn = .false.
540 END IF
541*
542 IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
543*
544* Scale by DMAX
545*
546 temp = abs( d( 1 ) )
547 DO 110 i = 2, rank
548 temp = max( temp, abs( d( i ) ) )
549 110 CONTINUE
550*
551 IF( temp.GT.zero ) THEN
552 alpha = dmax / temp
553 ELSE
554 info = 2
555 RETURN
556 END IF
557*
558 CALL dscal( rank, alpha, d, 1 )
559*
560 END IF
561*
562* 3) Generate Banded Matrix using Givens rotations.
563* Also the special case of UUB=LLB=0
564*
565* Compute Addressing constants to cover all
566* storage formats. Whether GE, SY, GB, or SB,
567* upper or lower triangle or both,
568* the (i,j)-th element is in
569* A( i - ISKEW*j + IOFFST, j )
570*
571 IF( ipack.GT.4 ) THEN
572 ilda = lda - 1
573 iskew = 1
574 IF( ipack.GT.5 ) THEN
575 ioffst = uub + 1
576 ELSE
577 ioffst = 1
578 END IF
579 ELSE
580 ilda = lda
581 iskew = 0
582 ioffst = 0
583 END IF
584*
585* IPACKG is the format that the matrix is generated in. If this is
586* different from IPACK, then the matrix must be repacked at the
587* end. It also signals how to compute the norm, for scaling.
588*
589 ipackg = 0
590 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
591*
592* Diagonal Matrix -- We are done, unless it
593* is to be stored SP/PP/TP (PACK='R' or 'C')
594*
595 IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
596 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
597 IF( ipack.LE.2 .OR. ipack.GE.5 )
598 $ ipackg = ipack
599*
600 ELSE IF( givens ) THEN
601*
602* Check whether to use Givens rotations,
603* Householder transformations, or nothing.
604*
605 IF( isym.EQ.1 ) THEN
606*
607* Non-symmetric -- A = U D V
608*
609 IF( ipack.GT.4 ) THEN
610 ipackg = ipack
611 ELSE
612 ipackg = 0
613 END IF
614*
615 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
616*
617 IF( topdwn ) THEN
618 jkl = 0
619 DO 140 jku = 1, uub
620*
621* Transform from bandwidth JKL, JKU-1 to JKL, JKU
622*
623* Last row actually rotated is M
624* Last column actually rotated is MIN( M+JKU, N )
625*
626 DO 130 jr = 1, min( m+jku, n ) + jkl - 1
627 extra = zero
628 angle = twopi*dlarnd( 1, iseed )
629 c = cos( angle )
630 s = sin( angle )
631 icol = max( 1, jr-jkl )
632 IF( jr.LT.m ) THEN
633 il = min( n, jr+jku ) + 1 - icol
634 CALL dlarot( .true., jr.GT.jkl, .false., il, c,
635 $ s, a( jr-iskew*icol+ioffst, icol ),
636 $ ilda, extra, dummy )
637 END IF
638*
639* Chase "EXTRA" back up
640*
641 ir = jr
642 ic = icol
643 DO 120 jch = jr - jkl, 1, -jkl - jku
644 IF( ir.LT.m ) THEN
645 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
646 $ ic+1 ), extra, c, s, dummy )
647 END IF
648 irow = max( 1, jch-jku )
649 il = ir + 2 - irow
650 temp = zero
651 iltemp = jch.GT.jku
652 CALL dlarot( .false., iltemp, .true., il, c, -s,
653 $ a( irow-iskew*ic+ioffst, ic ),
654 $ ilda, temp, extra )
655 IF( iltemp ) THEN
656 CALL dlartg( a( irow+1-iskew*( ic+1 )+ioffst,
657 $ ic+1 ), temp, c, s, dummy )
658 icol = max( 1, jch-jku-jkl )
659 il = ic + 2 - icol
660 extra = zero
661 CALL dlarot( .true., jch.GT.jku+jkl, .true.,
662 $ il, c, -s, a( irow-iskew*icol+
663 $ ioffst, icol ), ilda, extra,
664 $ temp )
665 ic = icol
666 ir = irow
667 END IF
668 120 CONTINUE
669 130 CONTINUE
670 140 CONTINUE
671*
672 jku = uub
673 DO 170 jkl = 1, llb
674*
675* Transform from bandwidth JKL-1, JKU to JKL, JKU
676*
677 DO 160 jc = 1, min( n+jkl, m ) + jku - 1
678 extra = zero
679 angle = twopi*dlarnd( 1, iseed )
680 c = cos( angle )
681 s = sin( angle )
682 irow = max( 1, jc-jku )
683 IF( jc.LT.n ) THEN
684 il = min( m, jc+jkl ) + 1 - irow
685 CALL dlarot( .false., jc.GT.jku, .false., il, c,
686 $ s, a( irow-iskew*jc+ioffst, jc ),
687 $ ilda, extra, dummy )
688 END IF
689*
690* Chase "EXTRA" back up
691*
692 ic = jc
693 ir = irow
694 DO 150 jch = jc - jku, 1, -jkl - jku
695 IF( ic.LT.n ) THEN
696 CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
697 $ ic+1 ), extra, c, s, dummy )
698 END IF
699 icol = max( 1, jch-jkl )
700 il = ic + 2 - icol
701 temp = zero
702 iltemp = jch.GT.jkl
703 CALL dlarot( .true., iltemp, .true., il, c, -s,
704 $ a( ir-iskew*icol+ioffst, icol ),
705 $ ilda, temp, extra )
706 IF( iltemp ) THEN
707 CALL dlartg( a( ir+1-iskew*( icol+1 )+ioffst,
708 $ icol+1 ), temp, c, s, dummy )
709 irow = max( 1, jch-jkl-jku )
710 il = ir + 2 - irow
711 extra = zero
712 CALL dlarot( .false., jch.GT.jkl+jku, .true.,
713 $ il, c, -s, a( irow-iskew*icol+
714 $ ioffst, icol ), ilda, extra,
715 $ temp )
716 ic = icol
717 ir = irow
718 END IF
719 150 CONTINUE
720 160 CONTINUE
721 170 CONTINUE
722*
723 ELSE
724*
725* Bottom-Up -- Start at the bottom right.
726*
727 jkl = 0
728 DO 200 jku = 1, uub
729*
730* Transform from bandwidth JKL, JKU-1 to JKL, JKU
731*
732* First row actually rotated is M
733* First column actually rotated is MIN( M+JKU, N )
734*
735 iendch = min( m, n+jkl ) - 1
736 DO 190 jc = min( m+jku, n ) - 1, 1 - jkl, -1
737 extra = zero
738 angle = twopi*dlarnd( 1, iseed )
739 c = cos( angle )
740 s = sin( angle )
741 irow = max( 1, jc-jku+1 )
742 IF( jc.GT.0 ) THEN
743 il = min( m, jc+jkl+1 ) + 1 - irow
744 CALL dlarot( .false., .false., jc+jkl.LT.m, il,
745 $ c, s, a( irow-iskew*jc+ioffst,
746 $ jc ), ilda, dummy, extra )
747 END IF
748*
749* Chase "EXTRA" back down
750*
751 ic = jc
752 DO 180 jch = jc + jkl, iendch, jkl + jku
753 ilextr = ic.GT.0
754 IF( ilextr ) THEN
755 CALL dlartg( a( jch-iskew*ic+ioffst, ic ),
756 $ extra, c, s, dummy )
757 END IF
758 ic = max( 1, ic )
759 icol = min( n-1, jch+jku )
760 iltemp = jch + jku.LT.n
761 temp = zero
762 CALL dlarot( .true., ilextr, iltemp, icol+2-ic,
763 $ c, s, a( jch-iskew*ic+ioffst, ic ),
764 $ ilda, extra, temp )
765 IF( iltemp ) THEN
766 CALL dlartg( a( jch-iskew*icol+ioffst,
767 $ icol ), temp, c, s, dummy )
768 il = min( iendch, jch+jkl+jku ) + 2 - jch
769 extra = zero
770 CALL dlarot( .false., .true.,
771 $ jch+jkl+jku.LE.iendch, il, c, s,
772 $ a( jch-iskew*icol+ioffst,
773 $ icol ), ilda, temp, extra )
774 ic = icol
775 END IF
776 180 CONTINUE
777 190 CONTINUE
778 200 CONTINUE
779*
780 jku = uub
781 DO 230 jkl = 1, llb
782*
783* Transform from bandwidth JKL-1, JKU to JKL, JKU
784*
785* First row actually rotated is MIN( N+JKL, M )
786* First column actually rotated is N
787*
788 iendch = min( n, m+jku ) - 1
789 DO 220 jr = min( n+jkl, m ) - 1, 1 - jku, -1
790 extra = zero
791 angle = twopi*dlarnd( 1, iseed )
792 c = cos( angle )
793 s = sin( angle )
794 icol = max( 1, jr-jkl+1 )
795 IF( jr.GT.0 ) THEN
796 il = min( n, jr+jku+1 ) + 1 - icol
797 CALL dlarot( .true., .false., jr+jku.LT.n, il,
798 $ c, s, a( jr-iskew*icol+ioffst,
799 $ icol ), ilda, dummy, extra )
800 END IF
801*
802* Chase "EXTRA" back down
803*
804 ir = jr
805 DO 210 jch = jr + jku, iendch, jkl + jku
806 ilextr = ir.GT.0
807 IF( ilextr ) THEN
808 CALL dlartg( a( ir-iskew*jch+ioffst, jch ),
809 $ extra, c, s, dummy )
810 END IF
811 ir = max( 1, ir )
812 irow = min( m-1, jch+jkl )
813 iltemp = jch + jkl.LT.m
814 temp = zero
815 CALL dlarot( .false., ilextr, iltemp, irow+2-ir,
816 $ c, s, a( ir-iskew*jch+ioffst,
817 $ jch ), ilda, extra, temp )
818 IF( iltemp ) THEN
819 CALL dlartg( a( irow-iskew*jch+ioffst, jch ),
820 $ temp, c, s, dummy )
821 il = min( iendch, jch+jkl+jku ) + 2 - jch
822 extra = zero
823 CALL dlarot( .true., .true.,
824 $ jch+jkl+jku.LE.iendch, il, c, s,
825 $ a( irow-iskew*jch+ioffst, jch ),
826 $ ilda, temp, extra )
827 ir = irow
828 END IF
829 210 CONTINUE
830 220 CONTINUE
831 230 CONTINUE
832 END IF
833*
834 ELSE
835*
836* Symmetric -- A = U D U'
837*
838 ipackg = ipack
839 ioffg = ioffst
840*
841 IF( topdwn ) THEN
842*
843* Top-Down -- Generate Upper triangle only
844*
845 IF( ipack.GE.5 ) THEN
846 ipackg = 6
847 ioffg = uub + 1
848 ELSE
849 ipackg = 1
850 END IF
851 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
852*
853 DO 260 k = 1, uub
854 DO 250 jc = 1, n - 1
855 irow = max( 1, jc-k )
856 il = min( jc+1, k+2 )
857 extra = zero
858 temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
859 angle = twopi*dlarnd( 1, iseed )
860 c = cos( angle )
861 s = sin( angle )
862 CALL dlarot( .false., jc.GT.k, .true., il, c, s,
863 $ a( irow-iskew*jc+ioffg, jc ), ilda,
864 $ extra, temp )
865 CALL dlarot( .true., .true., .false.,
866 $ min( k, n-jc )+1, c, s,
867 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
868 $ temp, dummy )
869*
870* Chase EXTRA back up the matrix
871*
872 icol = jc
873 DO 240 jch = jc - k, 1, -k
874 CALL dlartg( a( jch+1-iskew*( icol+1 )+ioffg,
875 $ icol+1 ), extra, c, s, dummy )
876 temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
877 CALL dlarot( .true., .true., .true., k+2, c, -s,
878 $ a( ( 1-iskew )*jch+ioffg, jch ),
879 $ ilda, temp, extra )
880 irow = max( 1, jch-k )
881 il = min( jch+1, k+2 )
882 extra = zero
883 CALL dlarot( .false., jch.GT.k, .true., il, c,
884 $ -s, a( irow-iskew*jch+ioffg, jch ),
885 $ ilda, extra, temp )
886 icol = jch
887 240 CONTINUE
888 250 CONTINUE
889 260 CONTINUE
890*
891* If we need lower triangle, copy from upper. Note that
892* the order of copying is chosen to work for 'q' -> 'b'
893*
894 IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
895 DO 280 jc = 1, n
896 irow = ioffst - iskew*jc
897 DO 270 jr = jc, min( n, jc+uub )
898 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
899 270 CONTINUE
900 280 CONTINUE
901 IF( ipack.EQ.5 ) THEN
902 DO 300 jc = n - uub + 1, n
903 DO 290 jr = n + 2 - jc, uub + 1
904 a( jr, jc ) = zero
905 290 CONTINUE
906 300 CONTINUE
907 END IF
908 IF( ipackg.EQ.6 ) THEN
909 ipackg = ipack
910 ELSE
911 ipackg = 0
912 END IF
913 END IF
914 ELSE
915*
916* Bottom-Up -- Generate Lower triangle only
917*
918 IF( ipack.GE.5 ) THEN
919 ipackg = 5
920 IF( ipack.EQ.6 )
921 $ ioffg = 1
922 ELSE
923 ipackg = 2
924 END IF
925 CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
926*
927 DO 330 k = 1, uub
928 DO 320 jc = n - 1, 1, -1
929 il = min( n+1-jc, k+2 )
930 extra = zero
931 temp = a( 1+( 1-iskew )*jc+ioffg, jc )
932 angle = twopi*dlarnd( 1, iseed )
933 c = cos( angle )
934 s = -sin( angle )
935 CALL dlarot( .false., .true., n-jc.GT.k, il, c, s,
936 $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
937 $ temp, extra )
938 icol = max( 1, jc-k+1 )
939 CALL dlarot( .true., .false., .true., jc+2-icol, c,
940 $ s, a( jc-iskew*icol+ioffg, icol ),
941 $ ilda, dummy, temp )
942*
943* Chase EXTRA back down the matrix
944*
945 icol = jc
946 DO 310 jch = jc + k, n - 1, k
947 CALL dlartg( a( jch-iskew*icol+ioffg, icol ),
948 $ extra, c, s, dummy )
949 temp = a( 1+( 1-iskew )*jch+ioffg, jch )
950 CALL dlarot( .true., .true., .true., k+2, c, s,
951 $ a( jch-iskew*icol+ioffg, icol ),
952 $ ilda, extra, temp )
953 il = min( n+1-jch, k+2 )
954 extra = zero
955 CALL dlarot( .false., .true., n-jch.GT.k, il, c,
956 $ s, a( ( 1-iskew )*jch+ioffg, jch ),
957 $ ilda, temp, extra )
958 icol = jch
959 310 CONTINUE
960 320 CONTINUE
961 330 CONTINUE
962*
963* If we need upper triangle, copy from lower. Note that
964* the order of copying is chosen to work for 'b' -> 'q'
965*
966 IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
967 DO 350 jc = n, 1, -1
968 irow = ioffst - iskew*jc
969 DO 340 jr = jc, max( 1, jc-uub ), -1
970 a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
971 340 CONTINUE
972 350 CONTINUE
973 IF( ipack.EQ.6 ) THEN
974 DO 370 jc = 1, uub
975 DO 360 jr = 1, uub + 1 - jc
976 a( jr, jc ) = zero
977 360 CONTINUE
978 370 CONTINUE
979 END IF
980 IF( ipackg.EQ.5 ) THEN
981 ipackg = ipack
982 ELSE
983 ipackg = 0
984 END IF
985 END IF
986 END IF
987 END IF
988*
989 ELSE
990*
991* 4) Generate Banded Matrix by first
992* Rotating by random Unitary matrices,
993* then reducing the bandwidth using Householder
994* transformations.
995*
996* Note: we should get here only if LDA .ge. N
997*
998 IF( isym.EQ.1 ) THEN
999*
1000* Non-symmetric -- A = U D V
1001*
1002 CALL dlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
1003 $ iinfo )
1004 ELSE
1005*
1006* Symmetric -- A = U D U'
1007*
1008 CALL dlagsy( m, llb, d, a, lda, iseed, work, iinfo )
1009*
1010 END IF
1011 IF( iinfo.NE.0 ) THEN
1012 info = 3
1013 RETURN
1014 END IF
1015 END IF
1016*
1017* 5) Pack the matrix
1018*
1019 IF( ipack.NE.ipackg ) THEN
1020 IF( ipack.EQ.1 ) THEN
1021*
1022* 'U' -- Upper triangular, not packed
1023*
1024 DO 390 j = 1, m
1025 DO 380 i = j + 1, m
1026 a( i, j ) = zero
1027 380 CONTINUE
1028 390 CONTINUE
1029*
1030 ELSE IF( ipack.EQ.2 ) THEN
1031*
1032* 'L' -- Lower triangular, not packed
1033*
1034 DO 410 j = 2, m
1035 DO 400 i = 1, j - 1
1036 a( i, j ) = zero
1037 400 CONTINUE
1038 410 CONTINUE
1039*
1040 ELSE IF( ipack.EQ.3 ) THEN
1041*
1042* 'C' -- Upper triangle packed Columnwise.
1043*
1044 icol = 1
1045 irow = 0
1046 DO 430 j = 1, m
1047 DO 420 i = 1, j
1048 irow = irow + 1
1049 IF( irow.GT.lda ) THEN
1050 irow = 1
1051 icol = icol + 1
1052 END IF
1053 a( irow, icol ) = a( i, j )
1054 420 CONTINUE
1055 430 CONTINUE
1056*
1057 ELSE IF( ipack.EQ.4 ) THEN
1058*
1059* 'R' -- Lower triangle packed Columnwise.
1060*
1061 icol = 1
1062 irow = 0
1063 DO 450 j = 1, m
1064 DO 440 i = j, m
1065 irow = irow + 1
1066 IF( irow.GT.lda ) THEN
1067 irow = 1
1068 icol = icol + 1
1069 END IF
1070 a( irow, icol ) = a( i, j )
1071 440 CONTINUE
1072 450 CONTINUE
1073*
1074 ELSE IF( ipack.GE.5 ) THEN
1075*
1076* 'B' -- The lower triangle is packed as a band matrix.
1077* 'Q' -- The upper triangle is packed as a band matrix.
1078* 'Z' -- The whole matrix is packed as a band matrix.
1079*
1080 IF( ipack.EQ.5 )
1081 $ uub = 0
1082 IF( ipack.EQ.6 )
1083 $ llb = 0
1084*
1085 DO 470 j = 1, uub
1086 DO 460 i = min( j+llb, m ), 1, -1
1087 a( i-j+uub+1, j ) = a( i, j )
1088 460 CONTINUE
1089 470 CONTINUE
1090*
1091 DO 490 j = uub + 2, n
1092 DO 480 i = j - uub, min( j+llb, m )
1093 a( i-j+uub+1, j ) = a( i, j )
1094 480 CONTINUE
1095 490 CONTINUE
1096 END IF
1097*
1098* If packed, zero out extraneous elements.
1099*
1100* Symmetric/Triangular Packed --
1101* zero out everything after A(IROW,ICOL)
1102*
1103 IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1104 DO 510 jc = icol, m
1105 DO 500 jr = irow + 1, lda
1106 a( jr, jc ) = zero
1107 500 CONTINUE
1108 irow = 0
1109 510 CONTINUE
1110*
1111 ELSE IF( ipack.GE.5 ) THEN
1112*
1113* Packed Band --
1114* 1st row is now in A( UUB+2-j, j), zero above it
1115* m-th row is now in A( M+UUB-j,j), zero below it
1116* last non-zero diagonal is now in A( UUB+LLB+1,j ),
1117* zero below it, too.
1118*
1119 ir1 = uub + llb + 2
1120 ir2 = uub + m + 2
1121 DO 540 jc = 1, n
1122 DO 520 jr = 1, uub + 1 - jc
1123 a( jr, jc ) = zero
1124 520 CONTINUE
1125 DO 530 jr = max( 1, min( ir1, ir2-jc ) ), lda
1126 a( jr, jc ) = zero
1127 530 CONTINUE
1128 540 CONTINUE
1129 END IF
1130 END IF
1131*
1132 RETURN
1133*
1134* End of DLATMT
1135*
1136 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dlagge(m, n, kl, ku, d, a, lda, iseed, work, info)
DLAGGE
Definition dlagge.f:113
subroutine dlagsy(n, k, d, a, lda, iseed, work, info)
DLAGSY
Definition dlagsy.f:101
subroutine dlarot(lrows, lleft, lright, nl, c, s, a, lda, xleft, xright)
DLAROT
Definition dlarot.f:226
subroutine dlatm7(mode, cond, irsign, idist, iseed, d, n, rank, info)
DLATM7
Definition dlatm7.f:122
subroutine dlatmt(m, n, dist, iseed, sym, d, mode, cond, dmax, rank, kl, ku, pack, a, lda, work, info)
DLATMT
Definition dlatmt.f:331
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79