/* 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_ */