/*
NOTE: This is generated code. Look in Misc/lapack_lite for information on
remaking this file.
*/
#include "f2c.h"
#ifdef HAVE_CONFIG
#include "config.h"
#else
extern doublereal dlamch_(char *);
#define EPSILON dlamch_("Epsilon")
#define SAFEMINIMUM dlamch_("Safe minimum")
#define PRECISION dlamch_("Precision")
#define BASE dlamch_("Base")
#endif
extern doublereal dlapy2_(doublereal *x, doublereal *y);
/*
f2c knows the exact rules for precedence, and so omits parentheses where not
strictly necessary. Since this is generated code, we don't really care if
it's readable, and we know what is written is correct. So don't warn about
them.
*/
#if defined(__GNUC__)
#pragma GCC diagnostic ignored "-Wparentheses"
#endif
/* Table of constant values */
static integer c__1 = 1;
static real c_b172 = 0.f;
static real c_b173 = 1.f;
static integer c__0 = 0;
integer ieeeck_(integer *ispec, real *zero, real *one)
{
/* System generated locals */
integer ret_val;
/* Local variables */
static real nan1, nan2, nan3, nan4, nan5, nan6, neginf, posinf, negzro,
newzro;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
June 2010
Purpose
=======
IEEECK is called from the ILAENV to verify that Infinity and
possibly NaN arithmetic is safe (i.e. will not trap).
Arguments
=========
ISPEC (input) INTEGER
Specifies whether to test just for inifinity arithmetic
or whether to test for infinity and NaN arithmetic.
= 0: Verify infinity arithmetic only.
= 1: Verify infinity and NaN arithmetic.
ZERO (input) REAL
Must contain the value 0.0
This is passed to prevent the compiler from optimizing
away this code.
ONE (input) REAL
Must contain the value 1.0
This is passed to prevent the compiler from optimizing
away this code.
RETURN VALUE: INTEGER
= 0: Arithmetic failed to produce the correct answers
= 1: Arithmetic produced the correct answers
*/
ret_val = 1;
posinf = *one / *zero;
if (posinf <= *one) {
ret_val = 0;
return ret_val;
}
neginf = -(*one) / *zero;
if (neginf >= *zero) {
ret_val = 0;
return ret_val;
}
negzro = *one / (neginf + *one);
if (negzro != *zero) {
ret_val = 0;
return ret_val;
}
neginf = *one / negzro;
if (neginf >= *zero) {
ret_val = 0;
return ret_val;
}
newzro = negzro + *zero;
if (newzro != *zero) {
ret_val = 0;
return ret_val;
}
posinf = *one / newzro;
if (posinf <= *one) {
ret_val = 0;
return ret_val;
}
neginf *= posinf;
if (neginf >= *zero) {
ret_val = 0;
return ret_val;
}
posinf *= posinf;
if (posinf <= *one) {
ret_val = 0;
return ret_val;
}
/* Return if we were only asked to check infinity arithmetic */
if (*ispec == 0) {
return ret_val;
}
nan1 = posinf + neginf;
nan2 = posinf / neginf;
nan3 = posinf / posinf;
nan4 = posinf * *zero;
nan5 = neginf * negzro;
nan6 = nan5 * *zero;
if (nan1 == nan1) {
ret_val = 0;
return ret_val;
}
if (nan2 == nan2) {
ret_val = 0;
return ret_val;
}
if (nan3 == nan3) {
ret_val = 0;
return ret_val;
}
if (nan4 == nan4) {
ret_val = 0;
return ret_val;
}
if (nan5 == nan5) {
ret_val = 0;
return ret_val;
}
if (nan6 == nan6) {
ret_val = 0;
return ret_val;
}
return ret_val;
} /* ieeeck_ */
integer ilaclc_(integer *m, integer *n, complex *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1, i__2;
/* Local variables */
static integer i__;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILACLC scans A for its last non-zero column.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) COMPLEX array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*n == 0) {
ret_val = *n;
} else /* if(complicated condition) */ {
i__1 = *n * a_dim1 + 1;
i__2 = *m + *n * a_dim1;
if (a[i__1].r != 0.f || a[i__1].i != 0.f || (a[i__2].r != 0.f || a[
i__2].i != 0.f)) {
ret_val = *n;
} else {
/* Now scan each column from the end, returning with the first non-zero. */
for (ret_val = *n; ret_val >= 1; --ret_val) {
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = i__ + ret_val * a_dim1;
if (a[i__2].r != 0.f || a[i__2].i != 0.f) {
return ret_val;
}
}
}
}
}
return ret_val;
} /* ilaclc_ */
integer ilaclr_(integer *m, integer *n, complex *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1, i__2;
/* Local variables */
static integer i__, j;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILACLR scans A for its last non-zero row.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) COMPLEX array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*m == 0) {
ret_val = *m;
} else /* if(complicated condition) */ {
i__1 = *m + a_dim1;
i__2 = *m + *n * a_dim1;
if (a[i__1].r != 0.f || a[i__1].i != 0.f || (a[i__2].r != 0.f || a[
i__2].i != 0.f)) {
ret_val = *m;
} else {
/* Scan up each column tracking the last zero row seen. */
ret_val = 0;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i__ = *m; i__ >= 1; --i__) {
i__2 = i__ + j * a_dim1;
if (a[i__2].r != 0.f || a[i__2].i != 0.f) {
goto L10;
}
}
L10:
ret_val = max(ret_val,i__);
}
}
}
return ret_val;
} /* ilaclr_ */
integer iladlc_(integer *m, integer *n, doublereal *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1;
/* Local variables */
static integer i__;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILADLC scans A for its last non-zero column.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*n == 0) {
ret_val = *n;
} else if (a[*n * a_dim1 + 1] != 0. || a[*m + *n * a_dim1] != 0.) {
ret_val = *n;
} else {
/* Now scan each column from the end, returning with the first non-zero. */
for (ret_val = *n; ret_val >= 1; --ret_val) {
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
if (a[i__ + ret_val * a_dim1] != 0.) {
return ret_val;
}
}
}
}
return ret_val;
} /* iladlc_ */
integer iladlr_(integer *m, integer *n, doublereal *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1;
/* Local variables */
static integer i__, j;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILADLR scans A for its last non-zero row.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*m == 0) {
ret_val = *m;
} else if (a[*m + a_dim1] != 0. || a[*m + *n * a_dim1] != 0.) {
ret_val = *m;
} else {
/* Scan up each column tracking the last zero row seen. */
ret_val = 0;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i__ = *m; i__ >= 1; --i__) {
if (a[i__ + j * a_dim1] != 0.) {
goto L10;
}
}
L10:
ret_val = max(ret_val,i__);
}
}
return ret_val;
} /* iladlr_ */
integer ilaenv_(integer *ispec, char *name__, char *opts, integer *n1,
integer *n2, integer *n3, integer *n4, ftnlen name_len, ftnlen
opts_len)
{
/* System generated locals */
integer ret_val;
/* Local variables */
static integer i__;
static char c1[1], c2[2], c3[3], c4[2];
static integer ic, nb, iz, nx;
static logical cname;
static integer nbmin;
static logical sname;
extern integer ieeeck_(integer *, real *, real *);
static char subnam[6];
extern integer iparmq_(integer *, char *, char *, integer *, integer *,
integer *, integer *, ftnlen, ftnlen);
/*
-- LAPACK auxiliary routine (version 3.2.1) --
-- April 2009 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILAENV is called from the LAPACK routines to choose problem-dependent
parameters for the local environment. See ISPEC for a description of
the parameters.
ILAENV returns an INTEGER
if ILAENV >= 0: ILAENV returns the value of the parameter specified by ISPEC
if ILAENV < 0: if ILAENV = -k, the k-th argument had an illegal value.
This version provides a set of parameters which should give good,
but not optimal, performance on many of the currently available
computers. Users are encouraged to modify this subroutine to set
the tuning parameters for their particular machine using the option
and problem size information in the arguments.
This routine will not function correctly if it is converted to all
lower case. Converting it to all upper case is allowed.
Arguments
=========
ISPEC (input) INTEGER
Specifies the parameter to be returned as the value of
ILAENV.
= 1: the optimal blocksize; if this value is 1, an unblocked
algorithm will give the best performance.
= 2: the minimum block size for which the block routine
should be used; if the usable block size is less than
this value, an unblocked routine should be used.
= 3: the crossover point (in a block routine, for N less
than this value, an unblocked routine should be used)
= 4: the number of shifts, used in the nonsymmetric
eigenvalue routines (DEPRECATED)
= 5: the minimum column dimension for blocking to be used;
rectangular blocks must have dimension at least k by m,
where k is given by ILAENV(2,...) and m by ILAENV(5,...)
= 6: the crossover point for the SVD (when reducing an m by n
matrix to bidiagonal form, if max(m,n)/min(m,n) exceeds
this value, a QR factorization is used first to reduce
the matrix to a triangular form.)
= 7: the number of processors
= 8: the crossover point for the multishift QR method
for nonsymmetric eigenvalue problems (DEPRECATED)
= 9: maximum size of the subproblems at the bottom of the
computation tree in the divide-and-conquer algorithm
(used by xGELSD and xGESDD)
=10: ieee NaN arithmetic can be trusted not to trap
=11: infinity arithmetic can be trusted not to trap
12 <= ISPEC <= 16:
xHSEQR or one of its subroutines,
see IPARMQ for detailed explanation
NAME (input) CHARACTER*(*)
The name of the calling subroutine, in either upper case or
lower case.
OPTS (input) CHARACTER*(*)
The character options to the subroutine NAME, concatenated
into a single character string. For example, UPLO = 'U',
TRANS = 'T', and DIAG = 'N' for a triangular routine would
be specified as OPTS = 'UTN'.
N1 (input) INTEGER
N2 (input) INTEGER
N3 (input) INTEGER
N4 (input) INTEGER
Problem dimensions for the subroutine NAME; these may not all
be required.
Further Details
===============
The following conventions have been used when calling ILAENV from the
LAPACK routines:
1) OPTS is a concatenation of all of the character options to
subroutine NAME, in the same order that they appear in the
argument list for NAME, even if they are not used in determining
the value of the parameter specified by ISPEC.
2) The problem dimensions N1, N2, N3, N4 are specified in the order
that they appear in the argument list for NAME. N1 is used
first, N2 second, and so on, and unused problem dimensions are
passed a value of -1.
3) The parameter value returned by ILAENV is checked for validity in
the calling subroutine. For example, ILAENV is used to retrieve
the optimal blocksize for STRTRI as follows:
NB = ILAENV( 1, 'STRTRI', UPLO // DIAG, N, -1, -1, -1 )
IF( NB.LE.1 ) NB = MAX( 1, N )
=====================================================================
*/
switch (*ispec) {
case 1: goto L10;
case 2: goto L10;
case 3: goto L10;
case 4: goto L80;
case 5: goto L90;
case 6: goto L100;
case 7: goto L110;
case 8: goto L120;
case 9: goto L130;
case 10: goto L140;
case 11: goto L150;
case 12: goto L160;
case 13: goto L160;
case 14: goto L160;
case 15: goto L160;
case 16: goto L160;
}
/* Invalid value for ISPEC */
ret_val = -1;
return ret_val;
L10:
/* Convert NAME to upper case if the first character is lower case. */
ret_val = 1;
s_copy(subnam, name__, (ftnlen)6, name_len);
ic = *(unsigned char *)subnam;
iz = 'Z';
if (iz == 90 || iz == 122) {
/* ASCII character set */
if (ic >= 97 && ic <= 122) {
*(unsigned char *)subnam = (char) (ic - 32);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 97 && ic <= 122) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
}
/* L20: */
}
}
} else if (iz == 233 || iz == 169) {
/* EBCDIC character set */
if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >= 162 &&
ic <= 169) {
*(unsigned char *)subnam = (char) (ic + 64);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 129 && ic <= 137 || ic >= 145 && ic <= 153 || ic >=
162 && ic <= 169) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic + 64);
}
/* L30: */
}
}
} else if (iz == 218 || iz == 250) {
/* Prime machines: ASCII+128 */
if (ic >= 225 && ic <= 250) {
*(unsigned char *)subnam = (char) (ic - 32);
for (i__ = 2; i__ <= 6; ++i__) {
ic = *(unsigned char *)&subnam[i__ - 1];
if (ic >= 225 && ic <= 250) {
*(unsigned char *)&subnam[i__ - 1] = (char) (ic - 32);
}
/* L40: */
}
}
}
*(unsigned char *)c1 = *(unsigned char *)subnam;
sname = *(unsigned char *)c1 == 'S' || *(unsigned char *)c1 == 'D';
cname = *(unsigned char *)c1 == 'C' || *(unsigned char *)c1 == 'Z';
if (! (cname || sname)) {
return ret_val;
}
s_copy(c2, subnam + 1, (ftnlen)2, (ftnlen)2);
s_copy(c3, subnam + 3, (ftnlen)3, (ftnlen)3);
s_copy(c4, c3 + 1, (ftnlen)2, (ftnlen)2);
switch (*ispec) {
case 1: goto L50;
case 2: goto L60;
case 3: goto L70;
}
L50:
/*
ISPEC = 1: block size
In these examples, separate code is provided for setting NB for
real and complex. We assume that NB will take the same value in
single or double precision.
*/
nb = 1;
if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
} else if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3,
"RQF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)
3, (ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3)
== 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 32;
} else {
nb = 32;
}
} else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "PO", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
} else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nb = 32;
} else if (sname && s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
nb = 64;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
nb = 64;
} else if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nb = 32;
} else if (s_cmp(c3, "GST", (ftnlen)3, (ftnlen)3) == 0) {
nb = 64;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nb = 32;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nb = 32;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nb = 32;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nb = 32;
}
}
} else if (s_cmp(c2, "GB", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
if (*n4 <= 64) {
nb = 1;
} else {
nb = 32;
}
} else {
if (*n4 <= 64) {
nb = 1;
} else {
nb = 32;
}
}
}
} else if (s_cmp(c2, "PB", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
if (*n2 <= 64) {
nb = 1;
} else {
nb = 32;
}
} else {
if (*n2 <= 64) {
nb = 1;
} else {
nb = 32;
}
}
}
} else if (s_cmp(c2, "TR", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (s_cmp(c2, "LA", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "UUM", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nb = 64;
} else {
nb = 64;
}
}
} else if (sname && s_cmp(c2, "ST", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "EBZ", (ftnlen)3, (ftnlen)3) == 0) {
nb = 1;
}
}
ret_val = nb;
return ret_val;
L60:
/* ISPEC = 2: minimum block size */
nbmin = 2;
if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
{
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
} else if (s_cmp(c3, "TRI", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nbmin = 2;
} else {
nbmin = 2;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRF", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nbmin = 8;
} else {
nbmin = 8;
}
} else if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nbmin = 2;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nbmin = 2;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nbmin = 2;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nbmin = 2;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nbmin = 2;
}
} else if (*(unsigned char *)c3 == 'M') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nbmin = 2;
}
}
}
ret_val = nbmin;
return ret_val;
L70:
/* ISPEC = 3: crossover point */
nx = 0;
if (s_cmp(c2, "GE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "QRF", (ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "RQF", (
ftnlen)3, (ftnlen)3) == 0 || s_cmp(c3, "LQF", (ftnlen)3, (
ftnlen)3) == 0 || s_cmp(c3, "QLF", (ftnlen)3, (ftnlen)3) == 0)
{
if (sname) {
nx = 128;
} else {
nx = 128;
}
} else if (s_cmp(c3, "HRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nx = 128;
} else {
nx = 128;
}
} else if (s_cmp(c3, "BRD", (ftnlen)3, (ftnlen)3) == 0) {
if (sname) {
nx = 128;
} else {
nx = 128;
}
}
} else if (s_cmp(c2, "SY", (ftnlen)2, (ftnlen)2) == 0) {
if (sname && s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nx = 32;
}
} else if (cname && s_cmp(c2, "HE", (ftnlen)2, (ftnlen)2) == 0) {
if (s_cmp(c3, "TRD", (ftnlen)3, (ftnlen)3) == 0) {
nx = 32;
}
} else if (sname && s_cmp(c2, "OR", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nx = 128;
}
}
} else if (cname && s_cmp(c2, "UN", (ftnlen)2, (ftnlen)2) == 0) {
if (*(unsigned char *)c3 == 'G') {
if (s_cmp(c4, "QR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "RQ",
(ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "LQ", (ftnlen)2, (
ftnlen)2) == 0 || s_cmp(c4, "QL", (ftnlen)2, (ftnlen)2) ==
0 || s_cmp(c4, "HR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(
c4, "TR", (ftnlen)2, (ftnlen)2) == 0 || s_cmp(c4, "BR", (
ftnlen)2, (ftnlen)2) == 0) {
nx = 128;
}
}
}
ret_val = nx;
return ret_val;
L80:
/* ISPEC = 4: number of shifts (used by xHSEQR) */
ret_val = 6;
return ret_val;
L90:
/* ISPEC = 5: minimum column dimension (not used) */
ret_val = 2;
return ret_val;
L100:
/* ISPEC = 6: crossover point for SVD (used by xGELSS and xGESVD) */
ret_val = (integer) ((real) min(*n1,*n2) * 1.6f);
return ret_val;
L110:
/* ISPEC = 7: number of processors (not used) */
ret_val = 1;
return ret_val;
L120:
/* ISPEC = 8: crossover point for multishift (used by xHSEQR) */
ret_val = 50;
return ret_val;
L130:
/*
ISPEC = 9: maximum size of the subproblems at the bottom of the
computation tree in the divide-and-conquer algorithm
(used by xGELSD and xGESDD)
*/
ret_val = 25;
return ret_val;
L140:
/*
ISPEC = 10: ieee NaN arithmetic can be trusted not to trap
ILAENV = 0
*/
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__1, &c_b172, &c_b173);
}
return ret_val;
L150:
/*
ISPEC = 11: infinity arithmetic can be trusted not to trap
ILAENV = 0
*/
ret_val = 1;
if (ret_val == 1) {
ret_val = ieeeck_(&c__0, &c_b172, &c_b173);
}
return ret_val;
L160:
/* 12 <= ISPEC <= 16: xHSEQR or one of its subroutines. */
ret_val = iparmq_(ispec, name__, opts, n1, n2, n3, n4, name_len, opts_len)
;
return ret_val;
/* End of ILAENV */
} /* ilaenv_ */
integer ilaslc_(integer *m, integer *n, real *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1;
/* Local variables */
static integer i__;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILASLC scans A for its last non-zero column.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) REAL array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*n == 0) {
ret_val = *n;
} else if (a[*n * a_dim1 + 1] != 0.f || a[*m + *n * a_dim1] != 0.f) {
ret_val = *n;
} else {
/* Now scan each column from the end, returning with the first non-zero. */
for (ret_val = *n; ret_val >= 1; --ret_val) {
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
if (a[i__ + ret_val * a_dim1] != 0.f) {
return ret_val;
}
}
}
}
return ret_val;
} /* ilaslc_ */
integer ilaslr_(integer *m, integer *n, real *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1;
/* Local variables */
static integer i__, j;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILASLR scans A for its last non-zero row.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) REAL array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*m == 0) {
ret_val = *m;
} else if (a[*m + a_dim1] != 0.f || a[*m + *n * a_dim1] != 0.f) {
ret_val = *m;
} else {
/* Scan up each column tracking the last zero row seen. */
ret_val = 0;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i__ = *m; i__ >= 1; --i__) {
if (a[i__ + j * a_dim1] != 0.f) {
goto L10;
}
}
L10:
ret_val = max(ret_val,i__);
}
}
return ret_val;
} /* ilaslr_ */
integer ilazlc_(integer *m, integer *n, doublecomplex *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1, i__2;
/* Local variables */
static integer i__;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILAZLC scans A for its last non-zero column.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) COMPLEX*16 array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*n == 0) {
ret_val = *n;
} else /* if(complicated condition) */ {
i__1 = *n * a_dim1 + 1;
i__2 = *m + *n * a_dim1;
if (a[i__1].r != 0. || a[i__1].i != 0. || (a[i__2].r != 0. || a[i__2]
.i != 0.)) {
ret_val = *n;
} else {
/* Now scan each column from the end, returning with the first non-zero. */
for (ret_val = *n; ret_val >= 1; --ret_val) {
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
i__2 = i__ + ret_val * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
return ret_val;
}
}
}
}
}
return ret_val;
} /* ilazlc_ */
integer ilazlr_(integer *m, integer *n, doublecomplex *a, integer *lda)
{
/* System generated locals */
integer a_dim1, a_offset, ret_val, i__1, i__2;
/* Local variables */
static integer i__, j;
/*
-- LAPACK auxiliary routine (version 3.2.2) --
-- June 2010 --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
Purpose
=======
ILAZLR scans A for its last non-zero row.
Arguments
=========
M (input) INTEGER
The number of rows of the matrix A.
N (input) INTEGER
The number of columns of the matrix A.
A (input) COMPLEX*16 array, dimension (LDA,N)
The m by n matrix A.
LDA (input) INTEGER
The leading dimension of the array A. LDA >= max(1,M).
=====================================================================
Quick test for the common case where one corner is non-zero.
*/
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
/* Function Body */
if (*m == 0) {
ret_val = *m;
} else /* if(complicated condition) */ {
i__1 = *m + a_dim1;
i__2 = *m + *n * a_dim1;
if (a[i__1].r != 0. || a[i__1].i != 0. || (a[i__2].r != 0. || a[i__2]
.i != 0.)) {
ret_val = *m;
} else {
/* Scan up each column tracking the last zero row seen. */
ret_val = 0;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
for (i__ = *m; i__ >= 1; --i__) {
i__2 = i__ + j * a_dim1;
if (a[i__2].r != 0. || a[i__2].i != 0.) {
goto L10;
}
}
L10:
ret_val = max(ret_val,i__);
}
}
}
return ret_val;
} /* ilazlr_ */
integer iparmq_(integer *ispec, char *name__, char *opts, integer *n, integer
*ilo, integer *ihi, integer *lwork, ftnlen name_len, ftnlen opts_len)
{
/* System generated locals */
integer ret_val, i__1, i__2;
real r__1;
/* Local variables */
static integer nh, ns;
/*
-- LAPACK auxiliary routine (version 3.2) --
-- LAPACK is a software package provided by Univ. of Tennessee, --
-- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
November 2006
Purpose
=======
This program sets problem and machine dependent parameters
useful for xHSEQR and its subroutines. It is called whenever
ILAENV is called with 12 <= ISPEC <= 16
Arguments
=========
ISPEC (input) integer scalar
ISPEC specifies which tunable parameter IPARMQ should
return.
ISPEC=12: (INMIN) Matrices of order nmin or less
are sent directly to xLAHQR, the implicit
double shift QR algorithm. NMIN must be
at least 11.
ISPEC=13: (INWIN) Size of the deflation window.
This is best set greater than or equal to
the number of simultaneous shifts NS.
Larger matrices benefit from larger deflation
windows.
ISPEC=14: (INIBL) Determines when to stop nibbling and
invest in an (expensive) multi-shift QR sweep.
If the aggressive early deflation subroutine
finds LD converged eigenvalues from an order
NW deflation window and LD.GT.(NW*NIBBLE)/100,
then the next QR sweep is skipped and early
deflation is applied immediately to the
remaining active diagonal block. Setting
IPARMQ(ISPEC=14) = 0 causes TTQRE to skip a
multi-shift QR sweep whenever early deflation
finds a converged eigenvalue. Setting
IPARMQ(ISPEC=14) greater than or equal to 100
prevents TTQRE from skipping a multi-shift
QR sweep.
ISPEC=15: (NSHFTS) The number of simultaneous shifts in
a multi-shift QR iteration.
ISPEC=16: (IACC22) IPARMQ is set to 0, 1 or 2 with the
following meanings.
0: During the multi-shift QR sweep,
xLAQR5 does not accumulate reflections and
does not use matrix-matrix multiply to
update the far-from-diagonal matrix
entries.
1: During the multi-shift QR sweep,
xLAQR5 and/or xLAQRaccumulates reflections and uses
matrix-matrix multiply to update the
far-from-diagonal matrix entries.
2: During the multi-shift QR sweep.
xLAQR5 accumulates reflections and takes
advantage of 2-by-2 block structure during
matrix-matrix multiplies.
(If xTRMM is slower than xGEMM, then
IPARMQ(ISPEC=16)=1 may be more efficient than
IPARMQ(ISPEC=16)=2 despite the greater level of
arithmetic work implied by the latter choice.)
NAME (input) character string
Name of the calling subroutine
OPTS (input) character string
This is a concatenation of the string arguments to
TTQRE.
N (input) integer scalar
N is the order of the Hessenberg matrix H.
ILO (input) INTEGER
IHI (input) INTEGER
It is assumed that H is already upper triangular
in rows and columns 1:ILO-1 and IHI+1:N.
LWORK (input) integer scalar
The amount of workspace available.
Further Details
===============
Little is known about how best to choose these parameters.
It is possible to use different values of the parameters
for each of CHSEQR, DHSEQR, SHSEQR and ZHSEQR.
It is probably best to choose different parameters for
different matrices and different parameters at different
times during the iteration, but this has not been
implemented --- yet.
The best choices of most of the parameters depend
in an ill-understood way on the relative execution
rate of xLAQR3 and xLAQR5 and on the nature of each
particular eigenvalue problem. Experiment may be the
only practical way to determine which choices are most
effective.
Following is a list of default values supplied by IPARMQ.
These defaults may be adjusted in order to attain better
performance in any particular computational environment.
IPARMQ(ISPEC=12) The xLAHQR vs xLAQR0 crossover point.
Default: 75. (Must be at least 11.)
IPARMQ(ISPEC=13) Recommended deflation window size.
This depends on ILO, IHI and NS, the
number of simultaneous shifts returned
by IPARMQ(ISPEC=15). The default for
(IHI-ILO+1).LE.500 is NS. The default
for (IHI-ILO+1).GT.500 is 3*NS/2.
IPARMQ(ISPEC=14) Nibble crossover point. Default: 14.
IPARMQ(ISPEC=15) Number of simultaneous shifts, NS.
a multi-shift QR iteration.
If IHI-ILO+1 is ...
greater than ...but less ... the
or equal to ... than default is
0 30 NS = 2+
30 60 NS = 4+
60 150 NS = 10
150 590 NS = **
590 3000 NS = 64
3000 6000 NS = 128
6000 infinity NS = 256
(+) By default matrices of this order are
passed to the implicit double shift routine
xLAHQR. See IPARMQ(ISPEC=12) above. These
values of NS are used only in case of a rare
xLAHQR failure.
(**) The asterisks (**) indicate an ad-hoc
function increasing from 10 to 64.
IPARMQ(ISPEC=16) Select structured matrix multiply.
(See ISPEC=16 above for details.)
Default: 3.
================================================================
*/
if (*ispec == 15 || *ispec == 13 || *ispec == 16) {
/* ==== Set the number simultaneous shifts ==== */
nh = *ihi - *ilo + 1;
ns = 2;
if (nh >= 30) {
ns = 4;
}
if (nh >= 60) {
ns = 10;
}
if (nh >= 150) {
/* Computing MAX */
r__1 = log((real) nh) / log(2.f);
i__1 = 10, i__2 = nh / i_nint(&r__1);
ns = max(i__1,i__2);
}
if (nh >= 590) {
ns = 64;
}
if (nh >= 3000) {
ns = 128;
}
if (nh >= 6000) {
ns = 256;
}
/* Computing MAX */
i__1 = 2, i__2 = ns - ns % 2;
ns = max(i__1,i__2);
}
if (*ispec == 12) {
/*
===== Matrices of order smaller than NMIN get sent
. to xLAHQR, the classic double shift algorithm.
. This must be at least 11. ====
*/
ret_val = 75;
} else if (*ispec == 14) {
/*
==== INIBL: skip a multi-shift qr iteration and
. whenever aggressive early deflation finds
. at least (NIBBLE*(window size)/100) deflations. ====
*/
ret_val = 14;
} else if (*ispec == 15) {
/* ==== NSHFTS: The number of simultaneous shifts ===== */
ret_val = ns;
} else if (*ispec == 13) {
/* ==== NW: deflation window size. ==== */
if (nh <= 500) {
ret_val = ns;
} else {
ret_val = ns * 3 / 2;
}
} else if (*ispec == 16) {
/*
==== IACC22: Whether to accumulate reflections
. before updating the far-from-diagonal elements
. and whether to use 2-by-2 block structure while
. doing it. A small amount of work could be saved
. by making this choice dependent also upon the
. NH=IHI-ILO+1.
*/
ret_val = 0;
if (ns >= 14) {
ret_val = 1;
}
if (ns >= 14) {
ret_val = 2;
}
} else {
/* ===== invalid value of ispec ===== */
ret_val = -1;
}
/* ==== End of IPARMQ ==== */
return ret_val;
} /* iparmq_ */