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