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