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