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