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