#include "f2c.h"
#include "blaswrap.h"
/* Table of constant values */
static doublecomplex c_b1 = {0.,0.};
static integer c__3 = 3;
/* Subroutine */ int zlatm4_(integer *itype, integer *n, integer *nz1,
integer *nz2, logical *rsign, doublereal *amagn, doublereal *rcond,
doublereal *triang, integer *idist, integer *iseed, doublecomplex *a,
integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
doublereal d__1, d__2;
doublecomplex z__1, z__2;
/* Builtin functions */
double pow_dd(doublereal *, doublereal *), log(doublereal), exp(
doublereal), z_abs(doublecomplex *);
/* Local variables */
integer i__, k, jc, jd, jr, kbeg, isdb, kend, isde, klen;
doublereal alpha;
doublecomplex ctemp;
extern doublereal dlaran_(integer *);
extern /* Double Complex */ VOID zlarnd_(doublecomplex *, integer *,
integer *);
extern /* Subroutine */ int zlaset_(char *, integer *, integer *,
doublecomplex *, doublecomplex *, doublecomplex *, integer *);
/* -- LAPACK auxiliary test routine (version 3.1) -- */
/* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/* November 2006 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* ZLATM4 generates basic square matrices, which may later be */
/* multiplied by others in order to produce test matrices. It is */
/* intended mainly to be used to test the generalized eigenvalue */
/* routines. */
/* It first generates the diagonal and (possibly) subdiagonal, */
/* according to the value of ITYPE, NZ1, NZ2, RSIGN, AMAGN, and RCOND. */
/* It then fills in the upper triangle with random numbers, if TRIANG is */
/* non-zero. */
/* Arguments */
/* ========= */
/* ITYPE (input) INTEGER */
/* The "type" of matrix on the diagonal and sub-diagonal. */
/* If ITYPE < 0, then type abs(ITYPE) is generated and then */
/* swapped end for end (A(I,J) := A'(N-J,N-I).) See also */
/* the description of AMAGN and RSIGN. */
/* Special types: */
/* = 0: the zero matrix. */
/* = 1: the identity. */
/* = 2: a transposed Jordan block. */
/* = 3: If N is odd, then a k+1 x k+1 transposed Jordan block */
/* followed by a k x k identity block, where k=(N-1)/2. */
/* If N is even, then k=(N-2)/2, and a zero diagonal entry */
/* is tacked onto the end. */
/* Diagonal types. The diagonal consists of NZ1 zeros, then */
/* k=N-NZ1-NZ2 nonzeros. The subdiagonal is zero. ITYPE */
/* specifies the nonzero diagonal entries as follows: */
/* = 4: 1, ..., k */
/* = 5: 1, RCOND, ..., RCOND */
/* = 6: 1, ..., 1, RCOND */
/* = 7: 1, a, a^2, ..., a^(k-1)=RCOND */
/* = 8: 1, 1-d, 1-2*d, ..., 1-(k-1)*d=RCOND */
/* = 9: random numbers chosen from (RCOND,1) */
/* = 10: random numbers with distribution IDIST (see ZLARND.) */
/* N (input) INTEGER */
/* The order of the matrix. */
/* NZ1 (input) INTEGER */
/* If abs(ITYPE) > 3, then the first NZ1 diagonal entries will */
/* be zero. */
/* NZ2 (input) INTEGER */
/* If abs(ITYPE) > 3, then the last NZ2 diagonal entries will */
/* be zero. */
/* RSIGN (input) LOGICAL */
/* = .TRUE.: The diagonal and subdiagonal entries will be */
/* multiplied by random numbers of magnitude 1. */
/* = .FALSE.: The diagonal and subdiagonal entries will be */
/* left as they are (usually non-negative real.) */
/* AMAGN (input) DOUBLE PRECISION */
/* The diagonal and subdiagonal entries will be multiplied by */
/* AMAGN. */
/* RCOND (input) DOUBLE PRECISION */
/* If abs(ITYPE) > 4, then the smallest diagonal entry will be */
/* RCOND. RCOND must be between 0 and 1. */
/* TRIANG (input) DOUBLE PRECISION */
/* The entries above the diagonal will be random numbers with */
/* magnitude bounded by TRIANG (i.e., random numbers multiplied */
/* by TRIANG.) */
/* IDIST (input) INTEGER */
/* On entry, DIST specifies the type of distribution to be used */
/* to generate a random matrix . */
/* = 1: real and imaginary parts each UNIFORM( 0, 1 ) */
/* = 2: real and imaginary parts each UNIFORM( -1, 1 ) */
/* = 3: real and imaginary parts each NORMAL( 0, 1 ) */
/* = 4: complex number uniform in DISK( 0, 1 ) */
/* ISEED (input/output) INTEGER array, dimension (4) */
/* On entry ISEED specifies the seed of the random number */
/* generator. The values of ISEED are changed on exit, and can */
/* be used in the next call to ZLATM4 to continue the same */
/* random number sequence. */
/* Note: ISEED(4) should be odd, for the random number generator */
/* used at present. */
/* A (output) COMPLEX*16 array, dimension (LDA, N) */
/* Array to be computed. */
/* LDA (input) INTEGER */
/* Leading dimension of A. Must be at least 1 and at least N. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Parameter adjustments */
--iseed;
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*n <= 0) {
return 0;
}
zlaset_("Full", n, n, &c_b1, &c_b1, &a[a_offset], lda);
/* Insure a correct ISEED */
if (iseed[4] % 2 != 1) {
++iseed[4];
}
/* Compute diagonal and subdiagonal according to ITYPE, NZ1, NZ2, */
/* and RCOND */
if (*itype != 0) {
if (abs(*itype) >= 4) {
/* Computing MAX */
/* Computing MIN */
i__3 = *n, i__4 = *nz1 + 1;
i__1 = 1, i__2 = min(i__3,i__4);
kbeg = max(i__1,i__2);
/* Computing MAX */
/* Computing MIN */
i__3 = *n, i__4 = *n - *nz2;
i__1 = kbeg, i__2 = min(i__3,i__4);
kend = max(i__1,i__2);
klen = kend + 1 - kbeg;
} else {
kbeg = 1;
kend = *n;
klen = *n;
}
isdb = 1;
isde = 0;
switch (abs(*itype)) {
case 1: goto L10;
case 2: goto L30;
case 3: goto L50;
case 4: goto L80;
case 5: goto L100;
case 6: goto L120;
case 7: goto L140;
case 8: goto L160;
case 9: goto L180;
case 10: goto L200;
}
/* abs(ITYPE) = 1: Identity */
L10:
i__1 = *n;
for (jd = 1; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
a[i__2].r = 1., a[i__2].i = 0.;
/* L20: */
}
goto L220;
/* abs(ITYPE) = 2: Transposed Jordan block */
L30:
i__1 = *n - 1;
for (jd = 1; jd <= i__1; ++jd) {
i__2 = jd + 1 + jd * a_dim1;
a[i__2].r = 1., a[i__2].i = 0.;
/* L40: */
}
isdb = 1;
isde = *n - 1;
goto L220;
/* abs(ITYPE) = 3: Transposed Jordan block, followed by the */
/* identity. */
L50:
k = (*n - 1) / 2;
i__1 = k;
for (jd = 1; jd <= i__1; ++jd) {
i__2 = jd + 1 + jd * a_dim1;
a[i__2].r = 1., a[i__2].i = 0.;
/* L60: */
}
isdb = 1;
isde = k;
i__1 = (k << 1) + 1;
for (jd = k + 2; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
a[i__2].r = 1., a[i__2].i = 0.;
/* L70: */
}
goto L220;
/* abs(ITYPE) = 4: 1,...,k */
L80:
i__1 = kend;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
i__3 = jd - *nz1;
z__1.r = (doublereal) i__3, z__1.i = 0.;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L90: */
}
goto L220;
/* abs(ITYPE) = 5: One large D value: */
L100:
i__1 = kend;
for (jd = kbeg + 1; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
z__1.r = *rcond, z__1.i = 0.;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L110: */
}
i__1 = kbeg + kbeg * a_dim1;
a[i__1].r = 1., a[i__1].i = 0.;
goto L220;
/* abs(ITYPE) = 6: One small D value: */
L120:
i__1 = kend - 1;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
a[i__2].r = 1., a[i__2].i = 0.;
/* L130: */
}
i__1 = kend + kend * a_dim1;
z__1.r = *rcond, z__1.i = 0.;
a[i__1].r = z__1.r, a[i__1].i = z__1.i;
goto L220;
/* abs(ITYPE) = 7: Exponentially distributed D values: */
L140:
i__1 = kbeg + kbeg * a_dim1;
a[i__1].r = 1., a[i__1].i = 0.;
if (klen > 1) {
d__1 = 1. / (doublereal) (klen - 1);
alpha = pow_dd(rcond, &d__1);
i__1 = klen;
for (i__ = 2; i__ <= i__1; ++i__) {
i__2 = *nz1 + i__ + (*nz1 + i__) * a_dim1;
d__2 = (doublereal) (i__ - 1);
d__1 = pow_dd(&alpha, &d__2);
z__1.r = d__1, z__1.i = 0.;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L150: */
}
}
goto L220;
/* abs(ITYPE) = 8: Arithmetically distributed D values: */
L160:
i__1 = kbeg + kbeg * a_dim1;
a[i__1].r = 1., a[i__1].i = 0.;
if (klen > 1) {
alpha = (1. - *rcond) / (doublereal) (klen - 1);
i__1 = klen;
for (i__ = 2; i__ <= i__1; ++i__) {
i__2 = *nz1 + i__ + (*nz1 + i__) * a_dim1;
d__1 = (doublereal) (klen - i__) * alpha + *rcond;
z__1.r = d__1, z__1.i = 0.;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L170: */
}
}
goto L220;
/* abs(ITYPE) = 9: Randomly distributed D values on ( RCOND, 1): */
L180:
alpha = log(*rcond);
i__1 = kend;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
d__1 = exp(alpha * dlaran_(&iseed[1]));
a[i__2].r = d__1, a[i__2].i = 0.;
/* L190: */
}
goto L220;
/* abs(ITYPE) = 10: Randomly distributed D values from DIST */
L200:
i__1 = kend;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
zlarnd_(&z__1, idist, &iseed[1]);
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
/* L210: */
}
L220:
/* Scale by AMAGN */
i__1 = kend;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
i__3 = jd + jd * a_dim1;
d__1 = *amagn * a[i__3].r;
a[i__2].r = d__1, a[i__2].i = 0.;
/* L230: */
}
i__1 = isde;
for (jd = isdb; jd <= i__1; ++jd) {
i__2 = jd + 1 + jd * a_dim1;
i__3 = jd + 1 + jd * a_dim1;
d__1 = *amagn * a[i__3].r;
a[i__2].r = d__1, a[i__2].i = 0.;
/* L240: */
}
/* If RSIGN = .TRUE., assign random signs to diagonal and */
/* subdiagonal */
if (*rsign) {
i__1 = kend;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
if (a[i__2].r != 0.) {
zlarnd_(&z__1, &c__3, &iseed[1]);
ctemp.r = z__1.r, ctemp.i = z__1.i;
d__1 = z_abs(&ctemp);
z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__2 = jd + jd * a_dim1;
i__3 = jd + jd * a_dim1;
d__1 = a[i__3].r;
z__1.r = d__1 * ctemp.r, z__1.i = d__1 * ctemp.i;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
}
/* L250: */
}
i__1 = isde;
for (jd = isdb; jd <= i__1; ++jd) {
i__2 = jd + 1 + jd * a_dim1;
if (a[i__2].r != 0.) {
zlarnd_(&z__1, &c__3, &iseed[1]);
ctemp.r = z__1.r, ctemp.i = z__1.i;
d__1 = z_abs(&ctemp);
z__1.r = ctemp.r / d__1, z__1.i = ctemp.i / d__1;
ctemp.r = z__1.r, ctemp.i = z__1.i;
i__2 = jd + 1 + jd * a_dim1;
i__3 = jd + 1 + jd * a_dim1;
d__1 = a[i__3].r;
z__1.r = d__1 * ctemp.r, z__1.i = d__1 * ctemp.i;
a[i__2].r = z__1.r, a[i__2].i = z__1.i;
}
/* L260: */
}
}
/* Reverse if ITYPE < 0 */
if (*itype < 0) {
i__1 = (kbeg + kend - 1) / 2;
for (jd = kbeg; jd <= i__1; ++jd) {
i__2 = jd + jd * a_dim1;
ctemp.r = a[i__2].r, ctemp.i = a[i__2].i;
i__2 = jd + jd * a_dim1;
i__3 = kbeg + kend - jd + (kbeg + kend - jd) * a_dim1;
a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
i__2 = kbeg + kend - jd + (kbeg + kend - jd) * a_dim1;
a[i__2].r = ctemp.r, a[i__2].i = ctemp.i;
/* L270: */
}
i__1 = (*n - 1) / 2;
for (jd = 1; jd <= i__1; ++jd) {
i__2 = jd + 1 + jd * a_dim1;
ctemp.r = a[i__2].r, ctemp.i = a[i__2].i;
i__2 = jd + 1 + jd * a_dim1;
i__3 = *n + 1 - jd + (*n - jd) * a_dim1;
a[i__2].r = a[i__3].r, a[i__2].i = a[i__3].i;
i__2 = *n + 1 - jd + (*n - jd) * a_dim1;
a[i__2].r = ctemp.r, a[i__2].i = ctemp.i;
/* L280: */
}
}
}
/* Fill in upper triangle */
if (*triang != 0.) {
i__1 = *n;
for (jc = 2; jc <= i__1; ++jc) {
i__2 = jc - 1;
for (jr = 1; jr <= i__2; ++jr) {
i__3 = jr + jc * a_dim1;
zlarnd_(&z__2, idist, &iseed[1]);
z__1.r = *triang * z__2.r, z__1.i = *triang * z__2.i;
a[i__3].r = z__1.r, a[i__3].i = z__1.i;
/* L290: */
}
/* L300: */
}
}
return 0;
/* End of ZLATM4 */
} /* zlatm4_ */