SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlatms()

subroutine zlatms ( integer  m,
integer  n,
character  dist,
integer, dimension( 4 )  iseed,
character  sym,
double precision, dimension( * )  d,
integer  mode,
double precision  cond,
double precision  dmax,
integer  kl,
integer  ku,
character  pack,
complex*16, dimension( lda, * )  a,
integer  lda,
complex*16, dimension( * )  work,
integer  info 
)

Definition at line 1 of file zlatms.f.

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