ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
zlatms.f
Go to the documentation of this file.
1  SUBROUTINE zlatms( M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX,
2  $ KL, KU, PACK, A, LDA, WORK, INFO )
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 *
1169  END
max
#define max(A, B)
Definition: pcgemr.c:180
zlaghe
subroutine zlaghe(N, K, D, A, LDA, ISEED, WORK, INFO)
Definition: zlaghe.f:2
zlarot
subroutine zlarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
Definition: zlarot.f:3
zlagsy
subroutine zlagsy(N, K, D, A, LDA, ISEED, WORK, INFO)
Definition: zlagsy.f:2
zlatms
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
Definition: zlatms.f:3
dlatm1
subroutine dlatm1(MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
Definition: dlatm1.f:2
zlagge
subroutine zlagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
Definition: zlagge.f:2
min
#define min(A, B)
Definition: pcgemr.c:181