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