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