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