1 /* -- translated by f2c (version 19940927).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11 static doublereal c_b23 = 0.;
12 static integer c__0 = 0;
13 static doublereal c_b39 = 1.;
14
dlatme_(integer * n,char * dist,integer * iseed,doublereal * d,integer * mode,doublereal * cond,doublereal * dmax__,char * ei,char * rsign,char * upper,char * sim,doublereal * ds,integer * modes,doublereal * conds,integer * kl,integer * ku,doublereal * anorm,doublereal * a,integer * lda,doublereal * work,integer * info)15 /* Subroutine */ int dlatme_(integer *n, char *dist, integer *iseed,
16 doublereal *d, integer *mode, doublereal *cond, doublereal *dmax__,
17 char *ei, char *rsign, char *upper, char *sim, doublereal *ds,
18 integer *modes, doublereal *conds, integer *kl, integer *ku,
19 doublereal *anorm, doublereal *a, integer *lda, doublereal *work,
20 integer *info)
21 {
22 /* System generated locals */
23 integer a_dim1, a_offset, i__1, i__2;
24 doublereal d__1, d__2, d__3;
25
26 /* Local variables */
27 static logical bads;
28 extern /* Subroutine */ int dger_(integer *, integer *, doublereal *,
29 doublereal *, integer *, doublereal *, integer *, doublereal *,
30 integer *);
31 static integer isim;
32 static doublereal temp;
33 static logical badei;
34 static integer i, j;
35 static doublereal alpha;
36 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
37 integer *);
38 extern logical lsame_(char *, char *);
39 extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
40 doublereal *, doublereal *, integer *, doublereal *, integer *,
41 doublereal *, doublereal *, integer *);
42 static integer iinfo;
43 static doublereal tempa[1];
44 static integer icols;
45 static logical useei;
46 static integer idist;
47 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
48 doublereal *, integer *);
49 static integer irows;
50 extern /* Subroutine */ int dlatm1_(integer *, doublereal *, integer *,
51 integer *, integer *, doublereal *, integer *, integer *);
52 static integer ic, jc;
53 extern doublereal dlange_(char *, integer *, integer *, doublereal *,
54 integer *, doublereal *);
55 static integer ir, jr;
56 extern /* Subroutine */ int dlarge_(integer *, doublereal *, integer *,
57 integer *, doublereal *, integer *), dlarfg_(integer *,
58 doublereal *, doublereal *, integer *, doublereal *);
59 extern doublereal dlaran_(integer *);
60 extern /* Subroutine */ int dlaset_(char *, integer *, integer *,
61 doublereal *, doublereal *, doublereal *, integer *),
62 xerbla_(char *, integer *), dlarnv_(integer *, integer *,
63 integer *, doublereal *);
64 static integer irsign, iupper;
65 static doublereal xnorms;
66 static integer jcr;
67 static doublereal tau;
68
69
70 /* -- LAPACK test routine (version 2.0) --
71 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
72 Courant Institute, Argonne National Lab, and Rice University
73 September 30, 1994
74
75
76 Purpose
77 =======
78
79 DLATME generates random non-symmetric square matrices with
80 specified eigenvalues for testing LAPACK programs.
81
82 DLATME operates by applying the following sequence of
83 operations:
84
85 1. Set the diagonal to D, where D may be input or
86 computed according to MODE, COND, DMAX, and RSIGN
87 as described below.
88
89 2. If complex conjugate pairs are desired (MODE=0 and EI(1)='R',
90 or MODE=5), certain pairs of adjacent elements of D are
91 interpreted as the real and complex parts of a complex
92 conjugate pair; A thus becomes block diagonal, with 1x1
93 and 2x2 blocks.
94
95 3. If UPPER='T', the upper triangle of A is set to random values
96 out of distribution DIST.
97
98 4. If SIM='T', A is multiplied on the left by a random matrix
99 X, whose singular values are specified by DS, MODES, and
100 CONDS, and on the right by X inverse.
101
102 5. If KL < N-1, the lower bandwidth is reduced to KL using
103 Householder transformations. If KU < N-1, the upper
104 bandwidth is reduced to KU.
105
106 6. If ANORM is not negative, the matrix is scaled to have
107 maximum-element-norm ANORM.
108
109 (Note: since the matrix cannot be reduced beyond Hessenberg form,
110
111 no packing options are available.)
112
113 Arguments
114 =========
115
116 N - INTEGER
117 The number of columns (or rows) of A. Not modified.
118
119 DIST - CHARACTER*1
120 On entry, DIST specifies the type of distribution to be used
121
122 to generate the random eigen-/singular values, and for the
123 upper triangle (see UPPER).
124 'U' => UNIFORM( 0, 1 ) ( 'U' for uniform )
125 'S' => UNIFORM( -1, 1 ) ( 'S' for symmetric )
126 'N' => NORMAL( 0, 1 ) ( 'N' for normal )
127 Not modified.
128
129 ISEED - INTEGER array, dimension ( 4 )
130 On entry ISEED specifies the seed of the random number
131 generator. They should lie between 0 and 4095 inclusive,
132 and ISEED(4) should be odd. The random number generator
133 uses a linear congruential sequence limited to small
134 integers, and so should produce machine independent
135 random numbers. The values of ISEED are changed on
136 exit, and can be used in the next call to DLATME
137 to continue the same random number sequence.
138 Changed on exit.
139
140 D - DOUBLE PRECISION array, dimension ( N )
141 This array is used to specify the eigenvalues of A. If
142 MODE=0, then D is assumed to contain the eigenvalues (but
143 see the description of EI), otherwise they will be
144 computed according to MODE, COND, DMAX, and RSIGN and
145 placed in D.
146 Modified if MODE is nonzero.
147
148 MODE - INTEGER
149 On entry this describes how the eigenvalues are to
150 be specified:
151 MODE = 0 means use D (with EI) as input
152 MODE = 1 sets D(1)=1 and D(2:N)=1.0/COND
153 MODE = 2 sets D(1:N-1)=1 and D(N)=1.0/COND
154 MODE = 3 sets D(I)=COND**(-(I-1)/(N-1))
155 MODE = 4 sets D(i)=1 - (i-1)/(N-1)*(1 - 1/COND)
156 MODE = 5 sets D to random numbers in the range
157 ( 1/COND , 1 ) such that their logarithms
158 are uniformly distributed. Each odd-even pair
159 of elements will be either used as two real
160 eigenvalues or as the real and imaginary part
161 of a complex conjugate pair of eigenvalues;
162 the choice of which is done is random, with
163 50-50 probability, for each pair.
164 MODE = 6 set D to random numbers from same distribution
165 as the rest of the matrix.
166 MODE < 0 has the same meaning as ABS(MODE), except that
167 the order of the elements of D is reversed.
168 Thus if MODE is between 1 and 4, D has entries ranging
169 from 1 to 1/COND, if between -1 and -4, D has entries
170 ranging from 1/COND to 1,
171 Not modified.
172
173 COND - DOUBLE PRECISION
174 On entry, this is used as described under MODE above.
175 If used, it must be >= 1. Not modified.
176
177 DMAX - DOUBLE PRECISION
178 If MODE is neither -6, 0 nor 6, the contents of D, as
179 computed according to MODE and COND, will be scaled by
180 DMAX / max(abs(D(i))). Note that DMAX need not be
181 positive: if DMAX is negative (or zero), D will be
182 scaled by a negative number (or zero).
183 Not modified.
184
185 EI - CHARACTER*1 array, dimension ( N )
186 If MODE is 0, and EI(1) is not ' ' (space character),
187 this array specifies which elements of D (on input) are
188 real eigenvalues and which are the real and imaginary parts
189
190 of a complex conjugate pair of eigenvalues. The elements
191 of EI may then only have the values 'R' and 'I'. If
192 EI(j)='R' and EI(j+1)='I', then the j-th eigenvalue is
193 CMPLX( D(j) , D(j+1) ), and the (j+1)-th is the complex
194 conjugate thereof. If EI(j)=EI(j+1)='R', then the j-th
195 eigenvalue is D(j) (i.e., real). EI(1) may not be 'I',
196 nor may two adjacent elements of EI both have the value 'I'.
197
198 If MODE is not 0, then EI is ignored. If MODE is 0 and
199 EI(1)=' ', then the eigenvalues will all be real.
200 Not modified.
201
202 RSIGN - CHARACTER*1
203 If MODE is not 0, 6, or -6, and RSIGN='T', then the
204 elements of D, as computed according to MODE and COND, will
205
206 be multiplied by a random sign (+1 or -1). If RSIGN='F',
207 they will not be. RSIGN may only have the values 'T' or
208 'F'.
209 Not modified.
210
211 UPPER - CHARACTER*1
212 If UPPER='T', then the elements of A above the diagonal
213 (and above the 2x2 diagonal blocks, if A has complex
214 eigenvalues) will be set to random numbers out of DIST.
215 If UPPER='F', they will not. UPPER may only have the
216 values 'T' or 'F'.
217 Not modified.
218
219 SIM - CHARACTER*1
220 If SIM='T', then A will be operated on by a "similarity
221 transform", i.e., multiplied on the left by a matrix X and
222 on the right by X inverse. X = U S V, where U and V are
223 random unitary matrices and S is a (diagonal) matrix of
224 singular values specified by DS, MODES, and CONDS. If
225 SIM='F', then A will not be transformed.
226 Not modified.
227
228 DS - DOUBLE PRECISION array, dimension ( N )
229 This array is used to specify the singular values of X,
230 in the same way that D specifies the eigenvalues of A.
231 If MODE=0, the DS contains the singular values, which
232 may not be zero.
233 Modified if MODE is nonzero.
234
235 MODES - INTEGER
236 CONDS - DOUBLE PRECISION
237 Same as MODE and COND, but for specifying the diagonal
238 of S. MODES=-6 and +6 are not allowed (since they would
239 result in randomly ill-conditioned eigenvalues.)
240
241 KL - INTEGER
242 This specifies the lower bandwidth of the matrix. KL=1
243 specifies upper Hessenberg form. If KL is at least N-1,
244 then A will have full lower bandwidth. KL must be at
245 least 1.
246 Not modified.
247
248 KU - INTEGER
249 This specifies the upper bandwidth of the matrix. KU=1
250 specifies lower Hessenberg form. If KU is at least N-1,
251 then A will have full upper bandwidth; if KU and KL
252 are both at least N-1, then A will be dense. Only one of
253 KU and KL may be less than N-1. KU must be at least 1.
254 Not modified.
255
256 ANORM - DOUBLE PRECISION
257 If ANORM is not negative, then A will be scaled by a non-
258 negative real number to make the maximum-element-norm of A
259 to be ANORM.
260 Not modified.
261
262 A - DOUBLE PRECISION array, dimension ( LDA, N )
263 On exit A is the desired test matrix.
264 Modified.
265
266 LDA - INTEGER
267 LDA specifies the first dimension of A as declared in the
268 calling program. LDA must be at least N.
269 Not modified.
270
271 WORK - DOUBLE PRECISION array, dimension ( 3*N )
272 Workspace.
273 Modified.
274
275 INFO - INTEGER
276 Error code. On exit, INFO will be set to one of the
277 following values:
278 0 => normal return
279 -1 => N negative
280 -2 => DIST illegal string
281 -5 => MODE not in range -6 to 6
282 -6 => COND less than 1.0, and MODE neither -6, 0 nor 6
283 -8 => EI(1) is not ' ' or 'R', EI(j) is not 'R' or 'I', or
284
285 two adjacent elements of EI are 'I'.
286 -9 => RSIGN is not 'T' or 'F'
287 -10 => UPPER is not 'T' or 'F'
288 -11 => SIM is not 'T' or 'F'
289 -12 => MODES=0 and DS has a zero singular value.
290 -13 => MODES is not in the range -5 to 5.
291 -14 => MODES is nonzero and CONDS is less than 1.
292 -15 => KL is less than 1.
293 -16 => KU is less than 1, or KL and KU are both less than
294 N-1.
295 -19 => LDA is less than N.
296 1 => Error return from DLATM1 (computing D)
297 2 => Cannot scale to DMAX (max. eigenvalue is 0)
298 3 => Error return from DLATM1 (computing DS)
299 4 => Error return from DLARGE
300 5 => Zero singular value from DLATM1.
301
302 =====================================================================
303
304
305
306 1) Decode and Test the input parameters.
307 Initialize flags & seed.
308
309 Parameter adjustments */
310 --iseed;
311 --d;
312 --ei;
313 --ds;
314 a_dim1 = *lda;
315 a_offset = a_dim1 + 1;
316 a -= a_offset;
317 --work;
318
319 /* Function Body */
320 *info = 0;
321
322 /* Quick return if possible */
323
324 if (*n == 0) {
325 return 0;
326 }
327
328 /* Decode DIST */
329
330 if (lsame_(dist, "U")) {
331 idist = 1;
332 } else if (lsame_(dist, "S")) {
333 idist = 2;
334 } else if (lsame_(dist, "N")) {
335 idist = 3;
336 } else {
337 idist = -1;
338 }
339
340 /* Check EI */
341
342 useei = TRUE_;
343 badei = FALSE_;
344 if (lsame_(ei + 1, " ") || *mode != 0) {
345 useei = FALSE_;
346 } else {
347 if (lsame_(ei + 1, "R")) {
348 i__1 = *n;
349 for (j = 2; j <= i__1; ++j) {
350 if (lsame_(ei + j, "I")) {
351 if (lsame_(ei + (j - 1), "I")) {
352 badei = TRUE_;
353 }
354 } else {
355 if (! lsame_(ei + j, "R")) {
356 badei = TRUE_;
357 }
358 }
359 /* L10: */
360 }
361 } else {
362 badei = TRUE_;
363 }
364 }
365
366 /* Decode RSIGN */
367
368 if (lsame_(rsign, "T")) {
369 irsign = 1;
370 } else if (lsame_(rsign, "F")) {
371 irsign = 0;
372 } else {
373 irsign = -1;
374 }
375
376 /* Decode UPPER */
377
378 if (lsame_(upper, "T")) {
379 iupper = 1;
380 } else if (lsame_(upper, "F")) {
381 iupper = 0;
382 } else {
383 iupper = -1;
384 }
385
386 /* Decode SIM */
387
388 if (lsame_(sim, "T")) {
389 isim = 1;
390 } else if (lsame_(sim, "F")) {
391 isim = 0;
392 } else {
393 isim = -1;
394 }
395
396 /* Check DS, if MODES=0 and ISIM=1 */
397
398 bads = FALSE_;
399 if (*modes == 0 && isim == 1) {
400 i__1 = *n;
401 for (j = 1; j <= i__1; ++j) {
402 if (ds[j] == 0.) {
403 bads = TRUE_;
404 }
405 /* L20: */
406 }
407 }
408
409 /* Set INFO if an error */
410
411 if (*n < 0) {
412 *info = -1;
413 } else if (idist == -1) {
414 *info = -2;
415 } else if (abs(*mode) > 6) {
416 *info = -5;
417 } else if (*mode != 0 && abs(*mode) != 6 && *cond < 1.) {
418 *info = -6;
419 } else if (badei) {
420 *info = -8;
421 } else if (irsign == -1) {
422 *info = -9;
423 } else if (iupper == -1) {
424 *info = -10;
425 } else if (isim == -1) {
426 *info = -11;
427 } else if (bads) {
428 *info = -12;
429 } else if (isim == 1 && abs(*modes) > 5) {
430 *info = -13;
431 } else if (isim == 1 && *modes != 0 && *conds < 1.) {
432 *info = -14;
433 } else if (*kl < 1) {
434 *info = -15;
435 } else if (*ku < 1 || *ku < *n - 1 && *kl < *n - 1) {
436 *info = -16;
437 } else if (*lda < max(1,*n)) {
438 *info = -19;
439 }
440
441 if (*info != 0) {
442 i__1 = -(*info);
443 xerbla_("DLATME", &i__1);
444 return 0;
445 }
446
447 /* Initialize random number generator */
448
449 for (i = 1; i <= 4; ++i) {
450 iseed[i] = (i__1 = iseed[i], abs(i__1)) % 4096;
451 /* L30: */
452 }
453
454 if (iseed[4] % 2 != 1) {
455 ++iseed[4];
456 }
457
458 /* 2) Set up diagonal of A
459
460 Compute D according to COND and MODE */
461
462 dlatm1_(mode, cond, &irsign, &idist, &iseed[1], &d[1], n, &iinfo);
463 if (iinfo != 0) {
464 *info = 1;
465 return 0;
466 }
467 if (*mode != 0 && abs(*mode) != 6) {
468
469 /* Scale by DMAX */
470
471 temp = abs(d[1]);
472 i__1 = *n;
473 for (i = 2; i <= i__1; ++i) {
474 /* Computing MAX */
475 d__2 = temp, d__3 = (d__1 = d[i], abs(d__1));
476 temp = max(d__2,d__3);
477 /* L40: */
478 }
479
480 if (temp > 0.) {
481 alpha = *dmax__ / temp;
482 } else if (*dmax__ != 0.) {
483 *info = 2;
484 return 0;
485 } else {
486 alpha = 0.;
487 }
488
489 dscal_(n, &alpha, &d[1], &c__1);
490
491 }
492
493 dlaset_("Full", n, n, &c_b23, &c_b23, &a[a_offset], lda);
494 i__1 = *lda + 1;
495 dcopy_(n, &d[1], &c__1, &a[a_offset], &i__1);
496
497 /* Set up complex conjugate pairs */
498
499 if (*mode == 0) {
500 if (useei) {
501 i__1 = *n;
502 for (j = 2; j <= i__1; ++j) {
503 if (lsame_(ei + j, "I")) {
504 a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
505 a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
506 a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
507 }
508 /* L50: */
509 }
510 }
511
512 } else if (abs(*mode) == 5) {
513
514 i__1 = *n;
515 for (j = 2; j <= i__1; j += 2) {
516 if (dlaran_(&iseed[1]) > .5) {
517 a[j - 1 + j * a_dim1] = a[j + j * a_dim1];
518 a[j + (j - 1) * a_dim1] = -a[j + j * a_dim1];
519 a[j + j * a_dim1] = a[j - 1 + (j - 1) * a_dim1];
520 }
521 /* L60: */
522 }
523 }
524
525 /* 3) If UPPER='T', set upper triangle of A to random numbers.
526 (but don't modify the corners of 2x2 blocks.) */
527
528 if (iupper != 0) {
529 i__1 = *n;
530 for (jc = 2; jc <= i__1; ++jc) {
531 if (a[jc - 1 + jc * a_dim1] != 0.) {
532 jr = jc - 2;
533 } else {
534 jr = jc - 1;
535 }
536 dlarnv_(&idist, &iseed[1], &jr, &a[jc * a_dim1 + 1]);
537 /* L70: */
538 }
539 }
540
541 /* 4) If SIM='T', apply similarity transformation.
542
543 -1
544 Transform is X A X , where X = U S V, thus
545
546 it is U S V A V' (1/S) U' */
547
548 if (isim != 0) {
549
550 /* Compute S (singular values of the eigenvector matrix)
551 according to CONDS and MODES */
552
553 dlatm1_(modes, conds, &c__0, &c__0, &iseed[1], &ds[1], n, &iinfo);
554 if (iinfo != 0) {
555 *info = 3;
556 return 0;
557 }
558
559 /* Multiply by V and V' */
560
561 dlarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
562 if (iinfo != 0) {
563 *info = 4;
564 return 0;
565 }
566
567 /* Multiply by S and (1/S) */
568
569 i__1 = *n;
570 for (j = 1; j <= i__1; ++j) {
571 dscal_(n, &ds[j], &a[j + a_dim1], lda);
572 if (ds[j] != 0.) {
573 d__1 = 1. / ds[j];
574 dscal_(n, &d__1, &a[j * a_dim1 + 1], &c__1);
575 } else {
576 *info = 5;
577 return 0;
578 }
579 /* L80: */
580 }
581
582 /* Multiply by U and U' */
583
584 dlarge_(n, &a[a_offset], lda, &iseed[1], &work[1], &iinfo);
585 if (iinfo != 0) {
586 *info = 4;
587 return 0;
588 }
589 }
590
591 /* 5) Reduce the bandwidth. */
592
593 if (*kl < *n - 1) {
594
595 /* Reduce bandwidth -- kill column */
596
597 i__1 = *n - 1;
598 for (jcr = *kl + 1; jcr <= i__1; ++jcr) {
599 ic = jcr - *kl;
600 irows = *n + 1 - jcr;
601 icols = *n + *kl - jcr;
602
603 dcopy_(&irows, &a[jcr + ic * a_dim1], &c__1, &work[1], &c__1);
604 xnorms = work[1];
605 dlarfg_(&irows, &xnorms, &work[2], &c__1, &tau);
606 work[1] = 1.;
607
608 dgemv_("T", &irows, &icols, &c_b39, &a[jcr + (ic + 1) * a_dim1],
609 lda, &work[1], &c__1, &c_b23, &work[irows + 1], &c__1)
610 ;
611 d__1 = -tau;
612 dger_(&irows, &icols, &d__1, &work[1], &c__1, &work[irows + 1], &
613 c__1, &a[jcr + (ic + 1) * a_dim1], lda);
614
615 dgemv_("N", n, &irows, &c_b39, &a[jcr * a_dim1 + 1], lda, &work[1]
616 , &c__1, &c_b23, &work[irows + 1], &c__1);
617 d__1 = -tau;
618 dger_(n, &irows, &d__1, &work[irows + 1], &c__1, &work[1], &c__1,
619 &a[jcr * a_dim1 + 1], lda);
620
621 a[jcr + ic * a_dim1] = xnorms;
622 i__2 = irows - 1;
623 dlaset_("Full", &i__2, &c__1, &c_b23, &c_b23, &a[jcr + 1 + ic *
624 a_dim1], lda);
625 /* L90: */
626 }
627 } else if (*ku < *n - 1) {
628
629 /* Reduce upper bandwidth -- kill a row at a time. */
630
631 i__1 = *n - 1;
632 for (jcr = *ku + 1; jcr <= i__1; ++jcr) {
633 ir = jcr - *ku;
634 irows = *n + *ku - jcr;
635 icols = *n + 1 - jcr;
636
637 dcopy_(&icols, &a[ir + jcr * a_dim1], lda, &work[1], &c__1);
638 xnorms = work[1];
639 dlarfg_(&icols, &xnorms, &work[2], &c__1, &tau);
640 work[1] = 1.;
641
642 dgemv_("N", &irows, &icols, &c_b39, &a[ir + 1 + jcr * a_dim1],
643 lda, &work[1], &c__1, &c_b23, &work[icols + 1], &c__1)
644 ;
645 d__1 = -tau;
646 dger_(&irows, &icols, &d__1, &work[icols + 1], &c__1, &work[1], &
647 c__1, &a[ir + 1 + jcr * a_dim1], lda);
648
649 dgemv_("C", &icols, n, &c_b39, &a[jcr + a_dim1], lda, &work[1], &
650 c__1, &c_b23, &work[icols + 1], &c__1);
651 d__1 = -tau;
652 dger_(&icols, n, &d__1, &work[1], &c__1, &work[icols + 1], &c__1,
653 &a[jcr + a_dim1], lda);
654
655 a[ir + jcr * a_dim1] = xnorms;
656 i__2 = icols - 1;
657 dlaset_("Full", &c__1, &i__2, &c_b23, &c_b23, &a[ir + (jcr + 1) *
658 a_dim1], lda);
659 /* L100: */
660 }
661 }
662
663 /* Scale the matrix to have norm ANORM */
664
665 if (*anorm >= 0.) {
666 temp = dlange_("M", n, n, &a[a_offset], lda, tempa);
667 if (temp > 0.) {
668 alpha = *anorm / temp;
669 i__1 = *n;
670 for (j = 1; j <= i__1; ++j) {
671 dscal_(n, &alpha, &a[j * a_dim1 + 1], &c__1);
672 /* L110: */
673 }
674 }
675 }
676
677 return 0;
678
679 /* End of DLATME */
680
681 } /* dlatme_ */
682
683