1 SUBROUTINE CLATME( N, DIST, ISEED, D, MODE, COND, DMAX, EI,
2 $ RSIGN,
3 $ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM,
4 $ A,
5 $ LDA, WORK, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * June 2010
10 *
11 * .. Scalar Arguments ..
12 CHARACTER DIST, EI, RSIGN, SIM, UPPER
13 INTEGER INFO, KL, KU, LDA, MODE, MODES, N
14 REAL ANORM, COND, CONDS
15 COMPLEX DMAX
16 * ..
17 * .. Array Arguments ..
18 INTEGER ISEED( 4 )
19 REAL DS( * )
20 COMPLEX A( LDA, * ), D( * ), WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * CLATME generates random non-symmetric square matrices with
27 * specified eigenvalues for testing LAPACK programs.
28 *
29 * CLATME operates by applying the following sequence of
30 * operations:
31 *
32 * 1. Set the diagonal to D, where D may be input or
33 * computed according to MODE, COND, DMAX, and RSIGN
34 * as described below.
35 *
36 * 2. If UPPER='T', the upper triangle of A is set to random values
37 * out of distribution DIST.
38 *
39 * 3. If SIM='T', A is multiplied on the left by a random matrix
40 * X, whose singular values are specified by DS, MODES, and
41 * CONDS, and on the right by X inverse.
42 *
43 * 4. If KL < N-1, the lower bandwidth is reduced to KL using
44 * Householder transformations. If KU < N-1, the upper
45 * bandwidth is reduced to KU.
46 *
47 * 5. If ANORM is not negative, the matrix is scaled to have
48 * maximum-element-norm ANORM.
49 *
50 * (Note: since the matrix cannot be reduced beyond Hessenberg form,
51 * no packing options are available.)
52 *
53 * Arguments
54 * =========
55 *
56 * N (input) INTEGER
57 * The number of columns (or rows) of A. Not modified.
58 *
59 * DIST (input) CHARACTER*1
60 * On entry, DIST specifies the type of distribution to be used
61 * to generate the random eigen-/singular values, and on the
62 * upper triangle (see UPPER).
63 * 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
64 * 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
65 * 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
66 * 'D' => uniform on the complex disc |z| < 1.
67 * Not modified.
68 *
69 * ISEED (input/output) INTEGER array, dimension ( 4 )
70 * On entry ISEED specifies the seed of the random number
71 * generator. They should lie between 0 and 4095 inclusive,
72 * and ISEED(4) should be odd. The random number generator
73 * uses a linear congruential sequence limited to small
74 * integers, and so should produce machine independent
75 * random numbers. The values of ISEED are changed on
76 * exit, and can be used in the next call to CLATME
77 * to continue the same random number sequence.
78 * Changed on exit.
79 *
80 * D (input/output) COMPLEX array, dimension ( N )
81 * This array is used to specify the eigenvalues of A. If
82 * MODE=0, then D is assumed to contain the eigenvalues
83 * otherwise they will be computed according to MODE, COND,
84 * DMAX, and RSIGN and placed in D.
85 * Modified if MODE is nonzero.
86 *
87 * MODE (input) INTEGER
88 * On entry this describes how the eigenvalues are to
89 * be specified:
90 * MODE = 0 means use D as input
91 * MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
92 * MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
93 * MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
94 * MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
95 * MODE = 5 sets D to random numbers in the range
96 * ( 1/COND , 1 ) such that their logarithms
97 * are uniformly distributed.
98 * MODE = 6 set D to random numbers from same distribution
99 * as the rest of the matrix.
100 * MODE < 0 has the same meaning as ABS(MODE), except that
101 * the order of the elements of D is reversed.
102 * Thus if MODE is between 1 and 4, D has entries ranging
103 * from 1 to 1/COND, if between -1 and -4, D has entries
104 * ranging from 1/COND to 1,
105 * Not modified.
106 *
107 * COND (input) REAL
108 * On entry, this is used as described under MODE above.
109 * If used, it must be >= 1. Not modified.
110 *
111 * DMAX (input) COMPLEX
112 * If MODE is neither -6, 0 nor 6, the contents of D, as
113 * computed according to MODE and COND, will be scaled by
114 * DMAX / max(abs(D(i))). Note that DMAX need not be
115 * positive or real: if DMAX is negative or complex (or zero),
116 * D will be scaled by a negative or complex number (or zero).
117 * If RSIGN='F' then the largest (absolute) eigenvalue will be
118 * equal to DMAX.
119 * Not modified.
120 *
121 * EI (input) CHARACTER*1 array, dimension ( N )
122 * (ignored)
123 * Not modified.
124 *
125 * RSIGN (input) CHARACTER*1
126 * If MODE is not 0, 6, or -6, and RSIGN='T', then the
127 * elements of D, as computed according to MODE and COND, will
128 * be multiplied by a random complex number from the unit
129 * circle |z| = 1. If RSIGN='F', they will not be. RSIGN may
130 * only have the values 'T' or 'F'.
131 * Not modified.
132 *
133 * UPPER (input) CHARACTER*1
134 * If UPPER='T', then the elements of A above the diagonal
135 * will be set to random numbers out of DIST. If UPPER='F',
136 * they will not. UPPER may only have the values 'T' or 'F'.
137 * Not modified.
138 *
139 * SIM (input) CHARACTER*1
140 * If SIM='T', then A will be operated on by a "similarity
141 * transform", i.e., multiplied on the left by a matrix X and
142 * on the right by X inverse. X = U S V, where U and V are
143 * random unitary matrices and S is a (diagonal) matrix of
144 * singular values specified by DS, MODES, and CONDS. If
145 * SIM='F', then A will not be transformed.
146 * Not modified.
147 *
148 * DS (input/output) REAL array, dimension ( N )
149 * This array is used to specify the singular values of X,
150 * in the same way that D specifies the eigenvalues of A.
151 * If MODE=0, the DS contains the singular values, which
152 * may not be zero.
153 * Modified if MODE is nonzero.
154 *
155 * MODES (input) INTEGER
156 *
157 * CONDS (input) REAL
158 * Similar to MODE and COND, but for specifying the diagonal
159 * of S. MODES=-6 and +6 are not allowed (since they would
160 * result in randomly ill-conditioned eigenvalues.)
161 *
162 * KL (input) INTEGER
163 * This specifies the lower bandwidth of the matrix. KL=1
164 * specifies upper Hessenberg form. If KL is at least N-1,
165 * then A will have full lower bandwidth.
166 * Not modified.
167 *
168 * KU (input) INTEGER
169 * This specifies the upper bandwidth of the matrix. KU=1
170 * specifies lower Hessenberg form. If KU is at least N-1,
171 * then A will have full upper bandwidth; if KU and KL
172 * are both at least N-1, then A will be dense. Only one of
173 * KU and KL may be less than N-1.
174 * Not modified.
175 *
176 * ANORM (input) REAL
177 * If ANORM is not negative, then A will be scaled by a non-
178 * negative real number to make the maximum-element-norm of A
179 * to be ANORM.
180 * Not modified.
181 *
182 * A (output) COMPLEX array, dimension ( LDA, N )
183 * On exit A is the desired test matrix.
184 * Modified.
185 *
186 * LDA (input) INTEGER
187 * LDA specifies the first dimension of A as declared in the
188 * calling program. LDA must be at least M.
189 * Not modified.
190 *
191 * WORK (workspace) COMPLEX array, dimension ( 3*N )
192 * Workspace.
193 * Modified.
194 *
195 * INFO (output) INTEGER
196 * Error code. On exit, INFO will be set to one of the
197 * following values:
198 * 0 => normal return
199 * -1 => N negative
200 * -2 => DIST illegal string
201 * -5 => MODE not in range -6 to 6
202 * -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
203 * -9 => RSIGN is not 'T' or 'F'
204 * -10 => UPPER is not 'T' or 'F'
205 * -11 => SIM is not 'T' or 'F'
206 * -12 => MODES=0 and DS has a zero singular value.
207 * -13 => MODES is not in the range -5 to 5.
208 * -14 => MODES is nonzero and CONDS is less than 1.
209 * -15 => KL is less than 1.
210 * -16 => KU is less than 1, or KL and KU are both less than
211 * N-1.
212 * -19 => LDA is less than M.
213 * 1 => Error return from CLATM1 (computing D)
214 * 2 => Cannot scale to DMAX (max. eigenvalue is 0)
215 * 3 => Error return from SLATM1 (computing DS)
216 * 4 => Error return from CLARGE
217 * 5 => Zero singular value from SLATM1.
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222 REAL ZERO
223 PARAMETER ( ZERO = 0.0E+0 )
224 REAL ONE
225 PARAMETER ( ONE = 1.0E+0 )
226 COMPLEX CZERO
227 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
228 COMPLEX CONE
229 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
230 * ..
231 * .. Local Scalars ..
232 LOGICAL BADS
233 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
234 $ ISIM, IUPPER, J, JC, JCR
235 REAL RALPHA, TEMP
236 COMPLEX ALPHA, TAU, XNORMS
237 * ..
238 * .. Local Arrays ..
239 REAL TEMPA( 1 )
240 * ..
241 * .. External Functions ..
242 LOGICAL LSAME
243 REAL CLANGE
244 COMPLEX CLARND
245 EXTERNAL LSAME, CLANGE, CLARND
246 * ..
247 * .. External Subroutines ..
248 EXTERNAL CCOPY, CGEMV, CGERC, CLACGV, CLARFG, CLARGE,
249 $ CLARNV, CLATM1, CLASET, CSCAL, CSSCAL, SLATM1,
250 $ XERBLA
251 * ..
252 * .. Intrinsic Functions ..
253 INTRINSIC ABS, CONJG, MAX, MOD
254 * ..
255 * .. Executable Statements ..
256 *
257 * 1) Decode and Test the input parameters.
258 * Initialize flags & seed.
259 *
260 INFO = 0
261 *
262 * Quick return if possible
263 *
264 IF( N.EQ.0 )
265 $ RETURN
266 *
267 * Decode DIST
268 *
269 IF( LSAME( DIST, 'U' ) ) THEN
270 IDIST = 1
271 ELSE IF( LSAME( DIST, 'S' ) ) THEN
272 IDIST = 2
273 ELSE IF( LSAME( DIST, 'N' ) ) THEN
274 IDIST = 3
275 ELSE IF( LSAME( DIST, 'D' ) ) THEN
276 IDIST = 4
277 ELSE
278 IDIST = -1
279 END IF
280 *
281 * Decode RSIGN
282 *
283 IF( LSAME( RSIGN, 'T' ) ) THEN
284 IRSIGN = 1
285 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN
286 IRSIGN = 0
287 ELSE
288 IRSIGN = -1
289 END IF
290 *
291 * Decode UPPER
292 *
293 IF( LSAME( UPPER, 'T' ) ) THEN
294 IUPPER = 1
295 ELSE IF( LSAME( UPPER, 'F' ) ) THEN
296 IUPPER = 0
297 ELSE
298 IUPPER = -1
299 END IF
300 *
301 * Decode SIM
302 *
303 IF( LSAME( SIM, 'T' ) ) THEN
304 ISIM = 1
305 ELSE IF( LSAME( SIM, 'F' ) ) THEN
306 ISIM = 0
307 ELSE
308 ISIM = -1
309 END IF
310 *
311 * Check DS, if MODES=0 and ISIM=1
312 *
313 BADS = .FALSE.
314 IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN
315 DO 10 J = 1, N
316 IF( DS( J ).EQ.ZERO )
317 $ BADS = .TRUE.
318 10 CONTINUE
319 END IF
320 *
321 * Set INFO if an error
322 *
323 IF( N.LT.0 ) THEN
324 INFO = -1
325 ELSE IF( IDIST.EQ.-1 ) THEN
326 INFO = -2
327 ELSE IF( ABS( MODE ).GT.6 ) THEN
328 INFO = -5
329 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
330 $ THEN
331 INFO = -6
332 ELSE IF( IRSIGN.EQ.-1 ) THEN
333 INFO = -9
334 ELSE IF( IUPPER.EQ.-1 ) THEN
335 INFO = -10
336 ELSE IF( ISIM.EQ.-1 ) THEN
337 INFO = -11
338 ELSE IF( BADS ) THEN
339 INFO = -12
340 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN
341 INFO = -13
342 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN
343 INFO = -14
344 ELSE IF( KL.LT.1 ) THEN
345 INFO = -15
346 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN
347 INFO = -16
348 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
349 INFO = -19
350 END IF
351 *
352 IF( INFO.NE.0 ) THEN
353 CALL XERBLA( 'CLATME', -INFO )
354 RETURN
355 END IF
356 *
357 * Initialize random number generator
358 *
359 DO 20 I = 1, 4
360 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
361 20 CONTINUE
362 *
363 IF( MOD( ISEED( 4 ), 2 ).NE.1 )
364 $ ISEED( 4 ) = ISEED( 4 ) + 1
365 *
366 * 2) Set up diagonal of A
367 *
368 * Compute D according to COND and MODE
369 *
370 CALL CLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO )
371 IF( IINFO.NE.0 ) THEN
372 INFO = 1
373 RETURN
374 END IF
375 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
376 *
377 * Scale by DMAX
378 *
379 TEMP = ABS( D( 1 ) )
380 DO 30 I = 2, N
381 TEMP = MAX( TEMP, ABS( D( I ) ) )
382 30 CONTINUE
383 *
384 IF( TEMP.GT.ZERO ) THEN
385 ALPHA = DMAX / TEMP
386 ELSE
387 INFO = 2
388 RETURN
389 END IF
390 *
391 CALL CSCAL( N, ALPHA, D, 1 )
392 *
393 END IF
394 *
395 CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
396 CALL CCOPY( N, D, 1, A, LDA+1 )
397 *
398 * 3) If UPPER='T', set upper triangle of A to random numbers.
399 *
400 IF( IUPPER.NE.0 ) THEN
401 DO 40 JC = 2, N
402 CALL CLARNV( IDIST, ISEED, JC-1, A( 1, JC ) )
403 40 CONTINUE
404 END IF
405 *
406 * 4) If SIM='T', apply similarity transformation.
407 *
408 * -1
409 * Transform is X A X , where X = U S V, thus
410 *
411 * it is U S V A V' (1/S) U'
412 *
413 IF( ISIM.NE.0 ) THEN
414 *
415 * Compute S (singular values of the eigenvector matrix)
416 * according to CONDS and MODES
417 *
418 CALL SLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO )
419 IF( IINFO.NE.0 ) THEN
420 INFO = 3
421 RETURN
422 END IF
423 *
424 * Multiply by V and V'
425 *
426 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO )
427 IF( IINFO.NE.0 ) THEN
428 INFO = 4
429 RETURN
430 END IF
431 *
432 * Multiply by S and (1/S)
433 *
434 DO 50 J = 1, N
435 CALL CSSCAL( N, DS( J ), A( J, 1 ), LDA )
436 IF( DS( J ).NE.ZERO ) THEN
437 CALL CSSCAL( N, ONE / DS( J ), A( 1, J ), 1 )
438 ELSE
439 INFO = 5
440 RETURN
441 END IF
442 50 CONTINUE
443 *
444 * Multiply by U and U'
445 *
446 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO )
447 IF( IINFO.NE.0 ) THEN
448 INFO = 4
449 RETURN
450 END IF
451 END IF
452 *
453 * 5) Reduce the bandwidth.
454 *
455 IF( KL.LT.N-1 ) THEN
456 *
457 * Reduce bandwidth -- kill column
458 *
459 DO 60 JCR = KL + 1, N - 1
460 IC = JCR - KL
461 IROWS = N + 1 - JCR
462 ICOLS = N + KL - JCR
463 *
464 CALL CCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 )
465 XNORMS = WORK( 1 )
466 CALL CLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU )
467 TAU = CONJG( TAU )
468 WORK( 1 ) = CONE
469 ALPHA = CLARND( 5, ISEED )
470 *
471 CALL CGEMV( 'C', IROWS, ICOLS, CONE, A( JCR, IC+1 ), LDA,
472 $ WORK, 1, CZERO, WORK( IROWS+1 ), 1 )
473 CALL CGERC( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1,
474 $ A( JCR, IC+1 ), LDA )
475 *
476 CALL CGEMV( 'N', N, IROWS, CONE, A( 1, JCR ), LDA, WORK, 1,
477 $ CZERO, WORK( IROWS+1 ), 1 )
478 CALL CGERC( N, IROWS, -CONJG( TAU ), WORK( IROWS+1 ), 1,
479 $ WORK, 1, A( 1, JCR ), LDA )
480 *
481 A( JCR, IC ) = XNORMS
482 CALL CLASET( 'Full', IROWS-1, 1, CZERO, CZERO,
483 $ A( JCR+1, IC ), LDA )
484 *
485 CALL CSCAL( ICOLS+1, ALPHA, A( JCR, IC ), LDA )
486 CALL CSCAL( N, CONJG( ALPHA ), A( 1, JCR ), 1 )
487 60 CONTINUE
488 ELSE IF( KU.LT.N-1 ) THEN
489 *
490 * Reduce upper bandwidth -- kill a row at a time.
491 *
492 DO 70 JCR = KU + 1, N - 1
493 IR = JCR - KU
494 IROWS = N + KU - JCR
495 ICOLS = N + 1 - JCR
496 *
497 CALL CCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 )
498 XNORMS = WORK( 1 )
499 CALL CLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU )
500 TAU = CONJG( TAU )
501 WORK( 1 ) = CONE
502 CALL CLACGV( ICOLS-1, WORK( 2 ), 1 )
503 ALPHA = CLARND( 5, ISEED )
504 *
505 CALL CGEMV( 'N', IROWS, ICOLS, CONE, A( IR+1, JCR ), LDA,
506 $ WORK, 1, CZERO, WORK( ICOLS+1 ), 1 )
507 CALL CGERC( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1,
508 $ A( IR+1, JCR ), LDA )
509 *
510 CALL CGEMV( 'C', ICOLS, N, CONE, A( JCR, 1 ), LDA, WORK, 1,
511 $ CZERO, WORK( ICOLS+1 ), 1 )
512 CALL CGERC( ICOLS, N, -CONJG( TAU ), WORK, 1,
513 $ WORK( ICOLS+1 ), 1, A( JCR, 1 ), LDA )
514 *
515 A( IR, JCR ) = XNORMS
516 CALL CLASET( 'Full', 1, ICOLS-1, CZERO, CZERO,
517 $ A( IR, JCR+1 ), LDA )
518 *
519 CALL CSCAL( IROWS+1, ALPHA, A( IR, JCR ), 1 )
520 CALL CSCAL( N, CONJG( ALPHA ), A( JCR, 1 ), LDA )
521 70 CONTINUE
522 END IF
523 *
524 * Scale the matrix to have norm ANORM
525 *
526 IF( ANORM.GE.ZERO ) THEN
527 TEMP = CLANGE( 'M', N, N, A, LDA, TEMPA )
528 IF( TEMP.GT.ZERO ) THEN
529 RALPHA = ANORM / TEMP
530 DO 80 J = 1, N
531 CALL CSSCAL( N, RALPHA, A( 1, J ), 1 )
532 80 CONTINUE
533 END IF
534 END IF
535 *
536 RETURN
537 *
538 * End of CLATME
539 *
540 END
2 $ RSIGN,
3 $ UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM,
4 $ A,
5 $ LDA, WORK, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * June 2010
10 *
11 * .. Scalar Arguments ..
12 CHARACTER DIST, EI, RSIGN, SIM, UPPER
13 INTEGER INFO, KL, KU, LDA, MODE, MODES, N
14 REAL ANORM, COND, CONDS
15 COMPLEX DMAX
16 * ..
17 * .. Array Arguments ..
18 INTEGER ISEED( 4 )
19 REAL DS( * )
20 COMPLEX A( LDA, * ), D( * ), WORK( * )
21 * ..
22 *
23 * Purpose
24 * =======
25 *
26 * CLATME generates random non-symmetric square matrices with
27 * specified eigenvalues for testing LAPACK programs.
28 *
29 * CLATME operates by applying the following sequence of
30 * operations:
31 *
32 * 1. Set the diagonal to D, where D may be input or
33 * computed according to MODE, COND, DMAX, and RSIGN
34 * as described below.
35 *
36 * 2. If UPPER='T', the upper triangle of A is set to random values
37 * out of distribution DIST.
38 *
39 * 3. If SIM='T', A is multiplied on the left by a random matrix
40 * X, whose singular values are specified by DS, MODES, and
41 * CONDS, and on the right by X inverse.
42 *
43 * 4. If KL < N-1, the lower bandwidth is reduced to KL using
44 * Householder transformations. If KU < N-1, the upper
45 * bandwidth is reduced to KU.
46 *
47 * 5. If ANORM is not negative, the matrix is scaled to have
48 * maximum-element-norm ANORM.
49 *
50 * (Note: since the matrix cannot be reduced beyond Hessenberg form,
51 * no packing options are available.)
52 *
53 * Arguments
54 * =========
55 *
56 * N (input) INTEGER
57 * The number of columns (or rows) of A. Not modified.
58 *
59 * DIST (input) CHARACTER*1
60 * On entry, DIST specifies the type of distribution to be used
61 * to generate the random eigen-/singular values, and on the
62 * upper triangle (see UPPER).
63 * 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
64 * 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
65 * 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
66 * 'D' => uniform on the complex disc |z| < 1.
67 * Not modified.
68 *
69 * ISEED (input/output) INTEGER array, dimension ( 4 )
70 * On entry ISEED specifies the seed of the random number
71 * generator. They should lie between 0 and 4095 inclusive,
72 * and ISEED(4) should be odd. The random number generator
73 * uses a linear congruential sequence limited to small
74 * integers, and so should produce machine independent
75 * random numbers. The values of ISEED are changed on
76 * exit, and can be used in the next call to CLATME
77 * to continue the same random number sequence.
78 * Changed on exit.
79 *
80 * D (input/output) COMPLEX array, dimension ( N )
81 * This array is used to specify the eigenvalues of A. If
82 * MODE=0, then D is assumed to contain the eigenvalues
83 * otherwise they will be computed according to MODE, COND,
84 * DMAX, and RSIGN and placed in D.
85 * Modified if MODE is nonzero.
86 *
87 * MODE (input) INTEGER
88 * On entry this describes how the eigenvalues are to
89 * be specified:
90 * MODE = 0 means use D as input
91 * MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
92 * MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
93 * MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
94 * MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
95 * MODE = 5 sets D to random numbers in the range
96 * ( 1/COND , 1 ) such that their logarithms
97 * are uniformly distributed.
98 * MODE = 6 set D to random numbers from same distribution
99 * as the rest of the matrix.
100 * MODE < 0 has the same meaning as ABS(MODE), except that
101 * the order of the elements of D is reversed.
102 * Thus if MODE is between 1 and 4, D has entries ranging
103 * from 1 to 1/COND, if between -1 and -4, D has entries
104 * ranging from 1/COND to 1,
105 * Not modified.
106 *
107 * COND (input) REAL
108 * On entry, this is used as described under MODE above.
109 * If used, it must be >= 1. Not modified.
110 *
111 * DMAX (input) COMPLEX
112 * If MODE is neither -6, 0 nor 6, the contents of D, as
113 * computed according to MODE and COND, will be scaled by
114 * DMAX / max(abs(D(i))). Note that DMAX need not be
115 * positive or real: if DMAX is negative or complex (or zero),
116 * D will be scaled by a negative or complex number (or zero).
117 * If RSIGN='F' then the largest (absolute) eigenvalue will be
118 * equal to DMAX.
119 * Not modified.
120 *
121 * EI (input) CHARACTER*1 array, dimension ( N )
122 * (ignored)
123 * Not modified.
124 *
125 * RSIGN (input) CHARACTER*1
126 * If MODE is not 0, 6, or -6, and RSIGN='T', then the
127 * elements of D, as computed according to MODE and COND, will
128 * be multiplied by a random complex number from the unit
129 * circle |z| = 1. If RSIGN='F', they will not be. RSIGN may
130 * only have the values 'T' or 'F'.
131 * Not modified.
132 *
133 * UPPER (input) CHARACTER*1
134 * If UPPER='T', then the elements of A above the diagonal
135 * will be set to random numbers out of DIST. If UPPER='F',
136 * they will not. UPPER may only have the values 'T' or 'F'.
137 * Not modified.
138 *
139 * SIM (input) CHARACTER*1
140 * If SIM='T', then A will be operated on by a "similarity
141 * transform", i.e., multiplied on the left by a matrix X and
142 * on the right by X inverse. X = U S V, where U and V are
143 * random unitary matrices and S is a (diagonal) matrix of
144 * singular values specified by DS, MODES, and CONDS. If
145 * SIM='F', then A will not be transformed.
146 * Not modified.
147 *
148 * DS (input/output) REAL array, dimension ( N )
149 * This array is used to specify the singular values of X,
150 * in the same way that D specifies the eigenvalues of A.
151 * If MODE=0, the DS contains the singular values, which
152 * may not be zero.
153 * Modified if MODE is nonzero.
154 *
155 * MODES (input) INTEGER
156 *
157 * CONDS (input) REAL
158 * Similar to MODE and COND, but for specifying the diagonal
159 * of S. MODES=-6 and +6 are not allowed (since they would
160 * result in randomly ill-conditioned eigenvalues.)
161 *
162 * KL (input) INTEGER
163 * This specifies the lower bandwidth of the matrix. KL=1
164 * specifies upper Hessenberg form. If KL is at least N-1,
165 * then A will have full lower bandwidth.
166 * Not modified.
167 *
168 * KU (input) INTEGER
169 * This specifies the upper bandwidth of the matrix. KU=1
170 * specifies lower Hessenberg form. If KU is at least N-1,
171 * then A will have full upper bandwidth; if KU and KL
172 * are both at least N-1, then A will be dense. Only one of
173 * KU and KL may be less than N-1.
174 * Not modified.
175 *
176 * ANORM (input) REAL
177 * If ANORM is not negative, then A will be scaled by a non-
178 * negative real number to make the maximum-element-norm of A
179 * to be ANORM.
180 * Not modified.
181 *
182 * A (output) COMPLEX array, dimension ( LDA, N )
183 * On exit A is the desired test matrix.
184 * Modified.
185 *
186 * LDA (input) INTEGER
187 * LDA specifies the first dimension of A as declared in the
188 * calling program. LDA must be at least M.
189 * Not modified.
190 *
191 * WORK (workspace) COMPLEX array, dimension ( 3*N )
192 * Workspace.
193 * Modified.
194 *
195 * INFO (output) INTEGER
196 * Error code. On exit, INFO will be set to one of the
197 * following values:
198 * 0 => normal return
199 * -1 => N negative
200 * -2 => DIST illegal string
201 * -5 => MODE not in range -6 to 6
202 * -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
203 * -9 => RSIGN is not 'T' or 'F'
204 * -10 => UPPER is not 'T' or 'F'
205 * -11 => SIM is not 'T' or 'F'
206 * -12 => MODES=0 and DS has a zero singular value.
207 * -13 => MODES is not in the range -5 to 5.
208 * -14 => MODES is nonzero and CONDS is less than 1.
209 * -15 => KL is less than 1.
210 * -16 => KU is less than 1, or KL and KU are both less than
211 * N-1.
212 * -19 => LDA is less than M.
213 * 1 => Error return from CLATM1 (computing D)
214 * 2 => Cannot scale to DMAX (max. eigenvalue is 0)
215 * 3 => Error return from SLATM1 (computing DS)
216 * 4 => Error return from CLARGE
217 * 5 => Zero singular value from SLATM1.
218 *
219 * =====================================================================
220 *
221 * .. Parameters ..
222 REAL ZERO
223 PARAMETER ( ZERO = 0.0E+0 )
224 REAL ONE
225 PARAMETER ( ONE = 1.0E+0 )
226 COMPLEX CZERO
227 PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) )
228 COMPLEX CONE
229 PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) )
230 * ..
231 * .. Local Scalars ..
232 LOGICAL BADS
233 INTEGER I, IC, ICOLS, IDIST, IINFO, IR, IROWS, IRSIGN,
234 $ ISIM, IUPPER, J, JC, JCR
235 REAL RALPHA, TEMP
236 COMPLEX ALPHA, TAU, XNORMS
237 * ..
238 * .. Local Arrays ..
239 REAL TEMPA( 1 )
240 * ..
241 * .. External Functions ..
242 LOGICAL LSAME
243 REAL CLANGE
244 COMPLEX CLARND
245 EXTERNAL LSAME, CLANGE, CLARND
246 * ..
247 * .. External Subroutines ..
248 EXTERNAL CCOPY, CGEMV, CGERC, CLACGV, CLARFG, CLARGE,
249 $ CLARNV, CLATM1, CLASET, CSCAL, CSSCAL, SLATM1,
250 $ XERBLA
251 * ..
252 * .. Intrinsic Functions ..
253 INTRINSIC ABS, CONJG, MAX, MOD
254 * ..
255 * .. Executable Statements ..
256 *
257 * 1) Decode and Test the input parameters.
258 * Initialize flags & seed.
259 *
260 INFO = 0
261 *
262 * Quick return if possible
263 *
264 IF( N.EQ.0 )
265 $ RETURN
266 *
267 * Decode DIST
268 *
269 IF( LSAME( DIST, 'U' ) ) THEN
270 IDIST = 1
271 ELSE IF( LSAME( DIST, 'S' ) ) THEN
272 IDIST = 2
273 ELSE IF( LSAME( DIST, 'N' ) ) THEN
274 IDIST = 3
275 ELSE IF( LSAME( DIST, 'D' ) ) THEN
276 IDIST = 4
277 ELSE
278 IDIST = -1
279 END IF
280 *
281 * Decode RSIGN
282 *
283 IF( LSAME( RSIGN, 'T' ) ) THEN
284 IRSIGN = 1
285 ELSE IF( LSAME( RSIGN, 'F' ) ) THEN
286 IRSIGN = 0
287 ELSE
288 IRSIGN = -1
289 END IF
290 *
291 * Decode UPPER
292 *
293 IF( LSAME( UPPER, 'T' ) ) THEN
294 IUPPER = 1
295 ELSE IF( LSAME( UPPER, 'F' ) ) THEN
296 IUPPER = 0
297 ELSE
298 IUPPER = -1
299 END IF
300 *
301 * Decode SIM
302 *
303 IF( LSAME( SIM, 'T' ) ) THEN
304 ISIM = 1
305 ELSE IF( LSAME( SIM, 'F' ) ) THEN
306 ISIM = 0
307 ELSE
308 ISIM = -1
309 END IF
310 *
311 * Check DS, if MODES=0 and ISIM=1
312 *
313 BADS = .FALSE.
314 IF( MODES.EQ.0 .AND. ISIM.EQ.1 ) THEN
315 DO 10 J = 1, N
316 IF( DS( J ).EQ.ZERO )
317 $ BADS = .TRUE.
318 10 CONTINUE
319 END IF
320 *
321 * Set INFO if an error
322 *
323 IF( N.LT.0 ) THEN
324 INFO = -1
325 ELSE IF( IDIST.EQ.-1 ) THEN
326 INFO = -2
327 ELSE IF( ABS( MODE ).GT.6 ) THEN
328 INFO = -5
329 ELSE IF( ( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) .AND. COND.LT.ONE )
330 $ THEN
331 INFO = -6
332 ELSE IF( IRSIGN.EQ.-1 ) THEN
333 INFO = -9
334 ELSE IF( IUPPER.EQ.-1 ) THEN
335 INFO = -10
336 ELSE IF( ISIM.EQ.-1 ) THEN
337 INFO = -11
338 ELSE IF( BADS ) THEN
339 INFO = -12
340 ELSE IF( ISIM.EQ.1 .AND. ABS( MODES ).GT.5 ) THEN
341 INFO = -13
342 ELSE IF( ISIM.EQ.1 .AND. MODES.NE.0 .AND. CONDS.LT.ONE ) THEN
343 INFO = -14
344 ELSE IF( KL.LT.1 ) THEN
345 INFO = -15
346 ELSE IF( KU.LT.1 .OR. ( KU.LT.N-1 .AND. KL.LT.N-1 ) ) THEN
347 INFO = -16
348 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
349 INFO = -19
350 END IF
351 *
352 IF( INFO.NE.0 ) THEN
353 CALL XERBLA( 'CLATME', -INFO )
354 RETURN
355 END IF
356 *
357 * Initialize random number generator
358 *
359 DO 20 I = 1, 4
360 ISEED( I ) = MOD( ABS( ISEED( I ) ), 4096 )
361 20 CONTINUE
362 *
363 IF( MOD( ISEED( 4 ), 2 ).NE.1 )
364 $ ISEED( 4 ) = ISEED( 4 ) + 1
365 *
366 * 2) Set up diagonal of A
367 *
368 * Compute D according to COND and MODE
369 *
370 CALL CLATM1( MODE, COND, IRSIGN, IDIST, ISEED, D, N, IINFO )
371 IF( IINFO.NE.0 ) THEN
372 INFO = 1
373 RETURN
374 END IF
375 IF( MODE.NE.0 .AND. ABS( MODE ).NE.6 ) THEN
376 *
377 * Scale by DMAX
378 *
379 TEMP = ABS( D( 1 ) )
380 DO 30 I = 2, N
381 TEMP = MAX( TEMP, ABS( D( I ) ) )
382 30 CONTINUE
383 *
384 IF( TEMP.GT.ZERO ) THEN
385 ALPHA = DMAX / TEMP
386 ELSE
387 INFO = 2
388 RETURN
389 END IF
390 *
391 CALL CSCAL( N, ALPHA, D, 1 )
392 *
393 END IF
394 *
395 CALL CLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
396 CALL CCOPY( N, D, 1, A, LDA+1 )
397 *
398 * 3) If UPPER='T', set upper triangle of A to random numbers.
399 *
400 IF( IUPPER.NE.0 ) THEN
401 DO 40 JC = 2, N
402 CALL CLARNV( IDIST, ISEED, JC-1, A( 1, JC ) )
403 40 CONTINUE
404 END IF
405 *
406 * 4) If SIM='T', apply similarity transformation.
407 *
408 * -1
409 * Transform is X A X , where X = U S V, thus
410 *
411 * it is U S V A V' (1/S) U'
412 *
413 IF( ISIM.NE.0 ) THEN
414 *
415 * Compute S (singular values of the eigenvector matrix)
416 * according to CONDS and MODES
417 *
418 CALL SLATM1( MODES, CONDS, 0, 0, ISEED, DS, N, IINFO )
419 IF( IINFO.NE.0 ) THEN
420 INFO = 3
421 RETURN
422 END IF
423 *
424 * Multiply by V and V'
425 *
426 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO )
427 IF( IINFO.NE.0 ) THEN
428 INFO = 4
429 RETURN
430 END IF
431 *
432 * Multiply by S and (1/S)
433 *
434 DO 50 J = 1, N
435 CALL CSSCAL( N, DS( J ), A( J, 1 ), LDA )
436 IF( DS( J ).NE.ZERO ) THEN
437 CALL CSSCAL( N, ONE / DS( J ), A( 1, J ), 1 )
438 ELSE
439 INFO = 5
440 RETURN
441 END IF
442 50 CONTINUE
443 *
444 * Multiply by U and U'
445 *
446 CALL CLARGE( N, A, LDA, ISEED, WORK, IINFO )
447 IF( IINFO.NE.0 ) THEN
448 INFO = 4
449 RETURN
450 END IF
451 END IF
452 *
453 * 5) Reduce the bandwidth.
454 *
455 IF( KL.LT.N-1 ) THEN
456 *
457 * Reduce bandwidth -- kill column
458 *
459 DO 60 JCR = KL + 1, N - 1
460 IC = JCR - KL
461 IROWS = N + 1 - JCR
462 ICOLS = N + KL - JCR
463 *
464 CALL CCOPY( IROWS, A( JCR, IC ), 1, WORK, 1 )
465 XNORMS = WORK( 1 )
466 CALL CLARFG( IROWS, XNORMS, WORK( 2 ), 1, TAU )
467 TAU = CONJG( TAU )
468 WORK( 1 ) = CONE
469 ALPHA = CLARND( 5, ISEED )
470 *
471 CALL CGEMV( 'C', IROWS, ICOLS, CONE, A( JCR, IC+1 ), LDA,
472 $ WORK, 1, CZERO, WORK( IROWS+1 ), 1 )
473 CALL CGERC( IROWS, ICOLS, -TAU, WORK, 1, WORK( IROWS+1 ), 1,
474 $ A( JCR, IC+1 ), LDA )
475 *
476 CALL CGEMV( 'N', N, IROWS, CONE, A( 1, JCR ), LDA, WORK, 1,
477 $ CZERO, WORK( IROWS+1 ), 1 )
478 CALL CGERC( N, IROWS, -CONJG( TAU ), WORK( IROWS+1 ), 1,
479 $ WORK, 1, A( 1, JCR ), LDA )
480 *
481 A( JCR, IC ) = XNORMS
482 CALL CLASET( 'Full', IROWS-1, 1, CZERO, CZERO,
483 $ A( JCR+1, IC ), LDA )
484 *
485 CALL CSCAL( ICOLS+1, ALPHA, A( JCR, IC ), LDA )
486 CALL CSCAL( N, CONJG( ALPHA ), A( 1, JCR ), 1 )
487 60 CONTINUE
488 ELSE IF( KU.LT.N-1 ) THEN
489 *
490 * Reduce upper bandwidth -- kill a row at a time.
491 *
492 DO 70 JCR = KU + 1, N - 1
493 IR = JCR - KU
494 IROWS = N + KU - JCR
495 ICOLS = N + 1 - JCR
496 *
497 CALL CCOPY( ICOLS, A( IR, JCR ), LDA, WORK, 1 )
498 XNORMS = WORK( 1 )
499 CALL CLARFG( ICOLS, XNORMS, WORK( 2 ), 1, TAU )
500 TAU = CONJG( TAU )
501 WORK( 1 ) = CONE
502 CALL CLACGV( ICOLS-1, WORK( 2 ), 1 )
503 ALPHA = CLARND( 5, ISEED )
504 *
505 CALL CGEMV( 'N', IROWS, ICOLS, CONE, A( IR+1, JCR ), LDA,
506 $ WORK, 1, CZERO, WORK( ICOLS+1 ), 1 )
507 CALL CGERC( IROWS, ICOLS, -TAU, WORK( ICOLS+1 ), 1, WORK, 1,
508 $ A( IR+1, JCR ), LDA )
509 *
510 CALL CGEMV( 'C', ICOLS, N, CONE, A( JCR, 1 ), LDA, WORK, 1,
511 $ CZERO, WORK( ICOLS+1 ), 1 )
512 CALL CGERC( ICOLS, N, -CONJG( TAU ), WORK, 1,
513 $ WORK( ICOLS+1 ), 1, A( JCR, 1 ), LDA )
514 *
515 A( IR, JCR ) = XNORMS
516 CALL CLASET( 'Full', 1, ICOLS-1, CZERO, CZERO,
517 $ A( IR, JCR+1 ), LDA )
518 *
519 CALL CSCAL( IROWS+1, ALPHA, A( IR, JCR ), 1 )
520 CALL CSCAL( N, CONJG( ALPHA ), A( JCR, 1 ), LDA )
521 70 CONTINUE
522 END IF
523 *
524 * Scale the matrix to have norm ANORM
525 *
526 IF( ANORM.GE.ZERO ) THEN
527 TEMP = CLANGE( 'M', N, N, A, LDA, TEMPA )
528 IF( TEMP.GT.ZERO ) THEN
529 RALPHA = ANORM / TEMP
530 DO 80 J = 1, N
531 CALL CSSCAL( N, RALPHA, A( 1, J ), 1 )
532 80 CONTINUE
533 END IF
534 END IF
535 *
536 RETURN
537 *
538 * End of CLATME
539 *
540 END