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