ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
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  irow = 1
360  icol = 1
361 *
362  IF( ipack.EQ.5 .OR. ipack.EQ.6 ) THEN
363  minlda = uub + 1
364  ELSE IF( ipack.EQ.7 ) THEN
365  minlda = llb + uub + 1
366  ELSE
367  minlda = m
368  END IF
369 *
370 * Use Givens rotation method if bandwidth small enough,
371 * or if LDA is too small to store the matrix unpacked.
372 *
373  givens = .false.
374  IF( isym.EQ.1 ) THEN
375  IF( dble( llb+uub ).LT.0.3d0*dble( max( 1, mr+nc ) ) )
376  $ givens = .true.
377  ELSE
378  IF( 2*llb.LT.m )
379  $ givens = .true.
380  END IF
381  IF( lda.LT.m .AND. lda.GE.minlda )
382  $ givens = .true.
383 *
384 * Set INFO if an error
385 *
386  IF( m.LT.0 ) THEN
387  info = -1
388  ELSE IF( m.NE.n .AND. isym.NE.1 ) THEN
389  info = -1
390  ELSE IF( n.LT.0 ) THEN
391  info = -2
392  ELSE IF( idist.EQ.-1 ) THEN
393  info = -3
394  ELSE IF( isym.EQ.-1 ) THEN
395  info = -5
396  ELSE IF( abs( mode ).GT.6 ) THEN
397  info = -7
398  ELSE IF( ( mode.NE.0 .AND. abs( mode ).NE.6 ) .AND. cond.LT.one )
399  $ THEN
400  info = -8
401  ELSE IF( kl.LT.0 ) THEN
402  info = -10
403  ELSE IF( ku.LT.0 .OR. ( isym.NE.1 .AND. kl.NE.ku ) ) THEN
404  info = -11
405  ELSE IF( ipack.EQ.-1 .OR. ( isympk.EQ.1 .AND. isym.EQ.1 ) .OR.
406  $ ( isympk.EQ.2 .AND. isym.EQ.1 .AND. kl.GT.0 ) .OR.
407  $ ( isympk.EQ.3 .AND. isym.EQ.1 .AND. ku.GT.0 ) .OR.
408  $ ( isympk.NE.0 .AND. m.NE.n ) ) THEN
409  info = -12
410  ELSE IF( lda.LT.max( 1, minlda ) ) THEN
411  info = -14
412  END IF
413 *
414  IF( info.NE.0 ) THEN
415  CALL xerbla( 'DLATMS', -info )
416  RETURN
417  END IF
418 *
419 * Initialize random number generator
420 *
421  DO 10 i = 1, 4
422  iseed( i ) = mod( abs( iseed( i ) ), 4096 )
423  10 CONTINUE
424 *
425  IF( mod( iseed( 4 ), 2 ).NE.1 )
426  $ iseed( 4 ) = iseed( 4 ) + 1
427 *
428 * 2) Set up D if indicated.
429 *
430 * Compute D according to COND and MODE
431 *
432  CALL dlatm1( mode, cond, irsign, idist, iseed, d, mnmin, iinfo )
433  IF( iinfo.NE.0 ) THEN
434  info = 1
435  RETURN
436  END IF
437 *
438 * Choose Top-Down if D is (apparently) increasing,
439 * Bottom-Up if D is (apparently) decreasing.
440 *
441  IF( abs( d( 1 ) ).LE.abs( d( mnmin ) ) ) THEN
442  topdwn = .true.
443  ELSE
444  topdwn = .false.
445  END IF
446 *
447  IF( mode.NE.0 .AND. abs( mode ).NE.6 ) THEN
448 *
449 * Scale by DMAX
450 *
451  temp = abs( d( 1 ) )
452  DO 20 i = 2, mnmin
453  temp = max( temp, abs( d( i ) ) )
454  20 CONTINUE
455 *
456  IF( temp.GT.zero ) THEN
457  alpha = dmax / temp
458  ELSE
459  info = 2
460  RETURN
461  END IF
462 *
463  CALL dscal( mnmin, alpha, d, 1 )
464 *
465  END IF
466 *
467 * 3) Generate Banded Matrix using Givens rotations.
468 * Also the special case of UUB=LLB=0
469 *
470 * Compute Addressing constants to cover all
471 * storage formats. Whether GE, SY, GB, or SB,
472 * upper or lower triangle or both,
473 * the (i,j)-th element is in
474 * A( i - ISKEW*j + IOFFST, j )
475 *
476  IF( ipack.GT.4 ) THEN
477  ilda = lda - 1
478  iskew = 1
479  IF( ipack.GT.5 ) THEN
480  ioffst = uub + 1
481  ELSE
482  ioffst = 1
483  END IF
484  ELSE
485  ilda = lda
486  iskew = 0
487  ioffst = 0
488  END IF
489 *
490 * IPACKG is the format that the matrix is generated in. If this is
491 * different from IPACK, then the matrix must be repacked at the
492 * end. It also signals how to compute the norm, for scaling.
493 *
494  ipackg = 0
495  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
496 *
497 * Diagonal Matrix -- We are done, unless it
498 * is to be stored SP/PP/TP (PACK='R' or 'C')
499 *
500  IF( llb.EQ.0 .AND. uub.EQ.0 ) THEN
501  CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
502  IF( ipack.LE.2 .OR. ipack.GE.5 )
503  $ ipackg = ipack
504 *
505  ELSE IF( givens ) THEN
506 *
507 * Check whether to use Givens rotations,
508 * Householder transformations, or nothing.
509 *
510  IF( isym.EQ.1 ) THEN
511 *
512 * Non-symmetric -- A = U D V
513 *
514  IF( ipack.GT.4 ) THEN
515  ipackg = ipack
516  ELSE
517  ipackg = 0
518  END IF
519 *
520  CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffst, 1 ), ilda+1 )
521 *
522  IF( topdwn ) THEN
523  jkl = 0
524  DO 50 jku = 1, uub
525 *
526 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
527 *
528 * Last row actually rotated is M
529 * Last column actually rotated is MIN( M+JKU, N )
530 *
531  DO 40 jr = 1, min( m+jku, n ) + jkl - 1
532  extra = zero
533  angle = twopi*dlarnd( 1, iseed )
534  c = cos( angle )
535  s = sin( angle )
536  icol = max( 1, jr-jkl )
537  IF( jr.LT.m ) THEN
538  il = min( n, jr+jku ) + 1 - icol
539  CALL dlarot( .true., jr.GT.jkl, .false., il, c,
540  $ s, a( jr-iskew*icol+ioffst, icol ),
541  $ ilda, extra, dummy )
542  END IF
543 *
544 * Chase "EXTRA" back up
545 *
546  ir = jr
547  ic = icol
548  DO 30 jch = jr - jkl, 1, -jkl - jku
549  IF( ir.LT.m ) THEN
550  CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
551  $ ic+1 ), extra, c, s, dummy )
552  END IF
553  irow = max( 1, jch-jku )
554  il = ir + 2 - irow
555  temp = zero
556  iltemp = jch.GT.jku
557  CALL dlarot( .false., iltemp, .true., il, c, -s,
558  $ a( irow-iskew*ic+ioffst, ic ),
559  $ ilda, temp, extra )
560  IF( iltemp ) THEN
561  CALL dlartg( a( irow+1-iskew*( ic+1 )+ioffst,
562  $ ic+1 ), temp, c, s, dummy )
563  icol = max( 1, jch-jku-jkl )
564  il = ic + 2 - icol
565  extra = zero
566  CALL dlarot( .true., jch.GT.jku+jkl, .true.,
567  $ il, c, -s, a( irow-iskew*icol+
568  $ ioffst, icol ), ilda, extra,
569  $ temp )
570  ic = icol
571  ir = irow
572  END IF
573  30 CONTINUE
574  40 CONTINUE
575  50 CONTINUE
576 *
577  jku = uub
578  DO 80 jkl = 1, llb
579 *
580 * Transform from bandwidth JKL-1, JKU to JKL, JKU
581 *
582  DO 70 jc = 1, min( n+jkl, m ) + jku - 1
583  extra = zero
584  angle = twopi*dlarnd( 1, iseed )
585  c = cos( angle )
586  s = sin( angle )
587  irow = max( 1, jc-jku )
588  IF( jc.LT.n ) THEN
589  il = min( m, jc+jkl ) + 1 - irow
590  CALL dlarot( .false., jc.GT.jku, .false., il, c,
591  $ s, a( irow-iskew*jc+ioffst, jc ),
592  $ ilda, extra, dummy )
593  END IF
594 *
595 * Chase "EXTRA" back up
596 *
597  ic = jc
598  ir = irow
599  DO 60 jch = jc - jku, 1, -jkl - jku
600  IF( ic.LT.n ) THEN
601  CALL dlartg( a( ir+1-iskew*( ic+1 )+ioffst,
602  $ ic+1 ), extra, c, s, dummy )
603  END IF
604  icol = max( 1, jch-jkl )
605  il = ic + 2 - icol
606  temp = zero
607  iltemp = jch.GT.jkl
608  CALL dlarot( .true., iltemp, .true., il, c, -s,
609  $ a( ir-iskew*icol+ioffst, icol ),
610  $ ilda, temp, extra )
611  IF( iltemp ) THEN
612  CALL dlartg( a( ir+1-iskew*( icol+1 )+ioffst,
613  $ icol+1 ), temp, c, s, dummy )
614  irow = max( 1, jch-jkl-jku )
615  il = ir + 2 - irow
616  extra = zero
617  CALL dlarot( .false., jch.GT.jkl+jku, .true.,
618  $ il, c, -s, a( irow-iskew*icol+
619  $ ioffst, icol ), ilda, extra,
620  $ temp )
621  ic = icol
622  ir = irow
623  END IF
624  60 CONTINUE
625  70 CONTINUE
626  80 CONTINUE
627 *
628  ELSE
629 *
630 * Bottom-Up -- Start at the bottom right.
631 *
632  jkl = 0
633  DO 110 jku = 1, uub
634 *
635 * Transform from bandwidth JKL, JKU-1 to JKL, JKU
636 *
637 * First row actually rotated is M
638 * First column actually rotated is MIN( M+JKU, N )
639 *
640  iendch = min( m, n+jkl ) - 1
641  DO 100 jc = min( m+jku, n ) - 1, 1 - jkl, -1
642  extra = zero
643  angle = twopi*dlarnd( 1, iseed )
644  c = cos( angle )
645  s = sin( angle )
646  irow = max( 1, jc-jku+1 )
647  IF( jc.GT.0 ) THEN
648  il = min( m, jc+jkl+1 ) + 1 - irow
649  CALL dlarot( .false., .false., jc+jkl.LT.m, il,
650  $ c, s, a( irow-iskew*jc+ioffst,
651  $ jc ), ilda, dummy, extra )
652  END IF
653 *
654 * Chase "EXTRA" back down
655 *
656  ic = jc
657  DO 90 jch = jc + jkl, iendch, jkl + jku
658  ilextr = ic.GT.0
659  IF( ilextr ) THEN
660  CALL dlartg( a( jch-iskew*ic+ioffst, ic ),
661  $ extra, c, s, dummy )
662  END IF
663  ic = max( 1, ic )
664  icol = min( n-1, jch+jku )
665  iltemp = jch + jku.LT.n
666  temp = zero
667  CALL dlarot( .true., ilextr, iltemp, icol+2-ic,
668  $ c, s, a( jch-iskew*ic+ioffst, ic ),
669  $ ilda, extra, temp )
670  IF( iltemp ) THEN
671  CALL dlartg( a( jch-iskew*icol+ioffst,
672  $ icol ), temp, c, s, dummy )
673  il = min( iendch, jch+jkl+jku ) + 2 - jch
674  extra = zero
675  CALL dlarot( .false., .true.,
676  $ jch+jkl+jku.LE.iendch, il, c, s,
677  $ a( jch-iskew*icol+ioffst,
678  $ icol ), ilda, temp, extra )
679  ic = icol
680  END IF
681  90 CONTINUE
682  100 CONTINUE
683  110 CONTINUE
684 *
685  jku = uub
686  DO 140 jkl = 1, llb
687 *
688 * Transform from bandwidth JKL-1, JKU to JKL, JKU
689 *
690 * First row actually rotated is MIN( N+JKL, M )
691 * First column actually rotated is N
692 *
693  iendch = min( n, m+jku ) - 1
694  DO 130 jr = min( n+jkl, m ) - 1, 1 - jku, -1
695  extra = zero
696  angle = twopi*dlarnd( 1, iseed )
697  c = cos( angle )
698  s = sin( angle )
699  icol = max( 1, jr-jkl+1 )
700  IF( jr.GT.0 ) THEN
701  il = min( n, jr+jku+1 ) + 1 - icol
702  CALL dlarot( .true., .false., jr+jku.LT.n, il,
703  $ c, s, a( jr-iskew*icol+ioffst,
704  $ icol ), ilda, dummy, extra )
705  END IF
706 *
707 * Chase "EXTRA" back down
708 *
709  ir = jr
710  DO 120 jch = jr + jku, iendch, jkl + jku
711  ilextr = ir.GT.0
712  IF( ilextr ) THEN
713  CALL dlartg( a( ir-iskew*jch+ioffst, jch ),
714  $ extra, c, s, dummy )
715  END IF
716  ir = max( 1, ir )
717  irow = min( m-1, jch+jkl )
718  iltemp = jch + jkl.LT.m
719  temp = zero
720  CALL dlarot( .false., ilextr, iltemp, irow+2-ir,
721  $ c, s, a( ir-iskew*jch+ioffst,
722  $ jch ), ilda, extra, temp )
723  IF( iltemp ) THEN
724  CALL dlartg( a( irow-iskew*jch+ioffst, jch ),
725  $ temp, c, s, dummy )
726  il = min( iendch, jch+jkl+jku ) + 2 - jch
727  extra = zero
728  CALL dlarot( .true., .true.,
729  $ jch+jkl+jku.LE.iendch, il, c, s,
730  $ a( irow-iskew*jch+ioffst, jch ),
731  $ ilda, temp, extra )
732  ir = irow
733  END IF
734  120 CONTINUE
735  130 CONTINUE
736  140 CONTINUE
737  END IF
738 *
739  ELSE
740 *
741 * Symmetric -- A = U D U'
742 *
743  ipackg = ipack
744  ioffg = ioffst
745 *
746  IF( topdwn ) THEN
747 *
748 * Top-Down -- Generate Upper triangle only
749 *
750  IF( ipack.GE.5 ) THEN
751  ipackg = 6
752  ioffg = uub + 1
753  ELSE
754  ipackg = 1
755  END IF
756  CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
757 *
758  DO 170 k = 1, uub
759  DO 160 jc = 1, n - 1
760  irow = max( 1, jc-k )
761  il = min( jc+1, k+2 )
762  extra = zero
763  temp = a( jc-iskew*( jc+1 )+ioffg, jc+1 )
764  angle = twopi*dlarnd( 1, iseed )
765  c = cos( angle )
766  s = sin( angle )
767  CALL dlarot( .false., jc.GT.k, .true., il, c, s,
768  $ a( irow-iskew*jc+ioffg, jc ), ilda,
769  $ extra, temp )
770  CALL dlarot( .true., .true., .false.,
771  $ min( k, n-jc )+1, c, s,
772  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
773  $ temp, dummy )
774 *
775 * Chase EXTRA back up the matrix
776 *
777  icol = jc
778  DO 150 jch = jc - k, 1, -k
779  CALL dlartg( a( jch+1-iskew*( icol+1 )+ioffg,
780  $ icol+1 ), extra, c, s, dummy )
781  temp = a( jch-iskew*( jch+1 )+ioffg, jch+1 )
782  CALL dlarot( .true., .true., .true., k+2, c, -s,
783  $ a( ( 1-iskew )*jch+ioffg, jch ),
784  $ ilda, temp, extra )
785  irow = max( 1, jch-k )
786  il = min( jch+1, k+2 )
787  extra = zero
788  CALL dlarot( .false., jch.GT.k, .true., il, c,
789  $ -s, a( irow-iskew*jch+ioffg, jch ),
790  $ ilda, extra, temp )
791  icol = jch
792  150 CONTINUE
793  160 CONTINUE
794  170 CONTINUE
795 *
796 * If we need lower triangle, copy from upper. Note that
797 * the order of copying is chosen to work for 'q' -> 'b'
798 *
799  IF( ipack.NE.ipackg .AND. ipack.NE.3 ) THEN
800  DO 190 jc = 1, n
801  irow = ioffst - iskew*jc
802  DO 180 jr = jc, min( n, jc+uub )
803  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
804  180 CONTINUE
805  190 CONTINUE
806  IF( ipack.EQ.5 ) THEN
807  DO 210 jc = n - uub + 1, n
808  DO 200 jr = n + 2 - jc, uub + 1
809  a( jr, jc ) = zero
810  200 CONTINUE
811  210 CONTINUE
812  END IF
813  IF( ipackg.EQ.6 ) THEN
814  ipackg = ipack
815  ELSE
816  ipackg = 0
817  END IF
818  END IF
819  ELSE
820 *
821 * Bottom-Up -- Generate Lower triangle only
822 *
823  IF( ipack.GE.5 ) THEN
824  ipackg = 5
825  IF( ipack.EQ.6 )
826  $ ioffg = 1
827  ELSE
828  ipackg = 2
829  END IF
830  CALL dcopy( mnmin, d, 1, a( 1-iskew+ioffg, 1 ), ilda+1 )
831 *
832  DO 240 k = 1, uub
833  DO 230 jc = n - 1, 1, -1
834  il = min( n+1-jc, k+2 )
835  extra = zero
836  temp = a( 1+( 1-iskew )*jc+ioffg, jc )
837  angle = twopi*dlarnd( 1, iseed )
838  c = cos( angle )
839  s = -sin( angle )
840  CALL dlarot( .false., .true., n-jc.GT.k, il, c, s,
841  $ a( ( 1-iskew )*jc+ioffg, jc ), ilda,
842  $ temp, extra )
843  icol = max( 1, jc-k+1 )
844  CALL dlarot( .true., .false., .true., jc+2-icol, c,
845  $ s, a( jc-iskew*icol+ioffg, icol ),
846  $ ilda, dummy, temp )
847 *
848 * Chase EXTRA back down the matrix
849 *
850  icol = jc
851  DO 220 jch = jc + k, n - 1, k
852  CALL dlartg( a( jch-iskew*icol+ioffg, icol ),
853  $ extra, c, s, dummy )
854  temp = a( 1+( 1-iskew )*jch+ioffg, jch )
855  CALL dlarot( .true., .true., .true., k+2, c, s,
856  $ a( jch-iskew*icol+ioffg, icol ),
857  $ ilda, extra, temp )
858  il = min( n+1-jch, k+2 )
859  extra = zero
860  CALL dlarot( .false., .true., n-jch.GT.k, il, c,
861  $ s, a( ( 1-iskew )*jch+ioffg, jch ),
862  $ ilda, temp, extra )
863  icol = jch
864  220 CONTINUE
865  230 CONTINUE
866  240 CONTINUE
867 *
868 * If we need upper triangle, copy from lower. Note that
869 * the order of copying is chosen to work for 'b' -> 'q'
870 *
871  IF( ipack.NE.ipackg .AND. ipack.NE.4 ) THEN
872  DO 260 jc = n, 1, -1
873  irow = ioffst - iskew*jc
874  DO 250 jr = jc, max( 1, jc-uub ), -1
875  a( jr+irow, jc ) = a( jc-iskew*jr+ioffg, jr )
876  250 CONTINUE
877  260 CONTINUE
878  IF( ipack.EQ.6 ) THEN
879  DO 280 jc = 1, uub
880  DO 270 jr = 1, uub + 1 - jc
881  a( jr, jc ) = zero
882  270 CONTINUE
883  280 CONTINUE
884  END IF
885  IF( ipackg.EQ.5 ) THEN
886  ipackg = ipack
887  ELSE
888  ipackg = 0
889  END IF
890  END IF
891  END IF
892  END IF
893 *
894  ELSE
895 *
896 * 4) Generate Banded Matrix by first
897 * Rotating by random Unitary matrices,
898 * then reducing the bandwidth using Householder
899 * transformations.
900 *
901 * Note: we should get here only if LDA .ge. N
902 *
903  IF( isym.EQ.1 ) THEN
904 *
905 * Non-symmetric -- A = U D V
906 *
907  CALL dlagge( mr, nc, llb, uub, d, a, lda, iseed, work,
908  $ iinfo )
909  ELSE
910 *
911 * Symmetric -- A = U D U'
912 *
913  CALL dlagsy( m, llb, d, a, lda, iseed, work, iinfo )
914 *
915  END IF
916  IF( iinfo.NE.0 ) THEN
917  info = 3
918  RETURN
919  END IF
920  END IF
921 *
922 * 5) Pack the matrix
923 *
924  IF( ipack.NE.ipackg ) THEN
925  IF( ipack.EQ.1 ) THEN
926 *
927 * 'U' -- Upper triangular, not packed
928 *
929  DO 300 j = 1, m
930  DO 290 i = j + 1, m
931  a( i, j ) = zero
932  290 CONTINUE
933  300 CONTINUE
934 *
935  ELSE IF( ipack.EQ.2 ) THEN
936 *
937 * 'L' -- Lower triangular, not packed
938 *
939  DO 320 j = 2, m
940  DO 310 i = 1, j - 1
941  a( i, j ) = zero
942  310 CONTINUE
943  320 CONTINUE
944 *
945  ELSE IF( ipack.EQ.3 ) THEN
946 *
947 * 'C' -- Upper triangle packed Columnwise.
948 *
949  icol = 1
950  irow = 0
951  DO 340 j = 1, m
952  DO 330 i = 1, j
953  irow = irow + 1
954  IF( irow.GT.lda ) THEN
955  irow = 1
956  icol = icol + 1
957  END IF
958  a( irow, icol ) = a( i, j )
959  330 CONTINUE
960  340 CONTINUE
961 *
962  ELSE IF( ipack.EQ.4 ) THEN
963 *
964 * 'R' -- Lower triangle packed Columnwise.
965 *
966  icol = 1
967  irow = 0
968  DO 360 j = 1, m
969  DO 350 i = j, m
970  irow = irow + 1
971  IF( irow.GT.lda ) THEN
972  irow = 1
973  icol = icol + 1
974  END IF
975  a( irow, icol ) = a( i, j )
976  350 CONTINUE
977  360 CONTINUE
978 *
979  ELSE IF( ipack.GE.5 ) THEN
980 *
981 * 'B' -- The lower triangle is packed as a band matrix.
982 * 'Q' -- The upper triangle is packed as a band matrix.
983 * 'Z' -- The whole matrix is packed as a band matrix.
984 *
985  IF( ipack.EQ.5 )
986  $ uub = 0
987  IF( ipack.EQ.6 )
988  $ llb = 0
989 *
990  DO 380 j = 1, uub
991  DO 370 i = min( j+llb, m ), 1, -1
992  a( i-j+uub+1, j ) = a( i, j )
993  370 CONTINUE
994  380 CONTINUE
995 *
996  DO 400 j = uub + 2, n
997  DO 390 i = j - uub, min( j+llb, m )
998  a( i-j+uub+1, j ) = a( i, j )
999  390 CONTINUE
1000  400 CONTINUE
1001  END IF
1002 *
1003 * If packed, zero out extraneous elements.
1004 *
1005 * Symmetric/Triangular Packed --
1006 * zero out everything after A(IROW,ICOL)
1007 *
1008  IF( ipack.EQ.3 .OR. ipack.EQ.4 ) THEN
1009  DO 420 jc = icol, m
1010  DO 410 jr = irow + 1, lda
1011  a( jr, jc ) = zero
1012  410 CONTINUE
1013  irow = 0
1014  420 CONTINUE
1015 *
1016  ELSE IF( ipack.GE.5 ) THEN
1017 *
1018 * Packed Band --
1019 * 1st row is now in A( UUB+2-j, j), zero above it
1020 * m-th row is now in A( M+UUB-j,j), zero below it
1021 * last non-zero diagonal is now in A( UUB+LLB+1,j ),
1022 * zero below it, too.
1023 *
1024  ir1 = uub + llb + 2
1025  ir2 = uub + m + 2
1026  DO 450 jc = 1, n
1027  DO 430 jr = 1, uub + 1 - jc
1028  a( jr, jc ) = zero
1029  430 CONTINUE
1030  DO 440 jr = max( 1, min( ir1, ir2-jc ) ), lda
1031  a( jr, jc ) = zero
1032  440 CONTINUE
1033  450 CONTINUE
1034  END IF
1035  END IF
1036 *
1037  RETURN
1038 *
1039 * End of DLATMS
1040 *
1041  END
max
#define max(A, B)
Definition: pcgemr.c:180
dlagge
subroutine dlagge(M, N, KL, KU, D, A, LDA, ISEED, WORK, INFO)
Definition: dlagge.f:2
dlarot
subroutine dlarot(LROWS, LLEFT, LRIGHT, NL, C, S, A, LDA, XLEFT, XRIGHT)
Definition: dlarot.f:3
dlagsy
subroutine dlagsy(N, K, D, A, LDA, ISEED, WORK, INFO)
Definition: dlagsy.f:2
dlatm1
subroutine dlatm1(MODE, COND, IRSIGN, IDIST, ISEED, D, N, INFO)
Definition: dlatm1.f:2
dlatms
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
Definition: dlatms.f:3
min
#define min(A, B)
Definition: pcgemr.c:181