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