From 35ea1c9256ee4b82b0fa68f1b51e30128547f5b9 Mon Sep 17 00:00:00 2001 From: Eric Wieser Date: Tue, 13 Dec 2016 19:41:52 +0000 Subject: MAINT: retranspile lapack from the supposed source The source used is http://archive.debian.org/debian/pool/main/l/lapack3/lapack3_3.0.20000531a.orig.tar.gz. Originally this was done with a patched f2c, but if the patch isn't provided in the source tree, there's no sensible way to use it --- numpy/linalg/lapack_lite/zlapack_lite.c | 4778 ++++++++++++------------------- 1 file changed, 1889 insertions(+), 2889 deletions(-) (limited to 'numpy/linalg/lapack_lite/zlapack_lite.c') diff --git a/numpy/linalg/lapack_lite/zlapack_lite.c b/numpy/linalg/lapack_lite/zlapack_lite.c index 29b017c89..0df4bdcda 100644 --- a/numpy/linalg/lapack_lite/zlapack_lite.c +++ b/numpy/linalg/lapack_lite/zlapack_lite.c @@ -233,9 +233,9 @@ L20: *info = -2; } else if (*n < 0) { *info = -3; - } else if ((*ilo < 1) || (*ilo > max(1,*n))) { + } else if (*ilo < 1 || *ilo > max(1,*n)) { *info = -4; - } else if ((*ihi < min(*ilo,*n)) || (*ihi > *n)) { + } else if (*ihi < min(*ilo,*n) || *ihi > *n) { *info = -5; } else if (*m < 0) { *info = -7; @@ -266,7 +266,7 @@ L20: /* Backward balance */ - if ((lsame_(job, "S")) || (lsame_(job, "B"))) { + if (lsame_(job, "S") || lsame_(job, "B")) { if (rightv) { i__1 = *ihi; @@ -296,7 +296,7 @@ L20: */ L30: - if ((lsame_(job, "P")) || (lsame_(job, "B"))) { + if (lsame_(job, "P") || lsame_(job, "B")) { if (rightv) { i__1 = *n; for (ii = 1; ii <= i__1; ++ii) { @@ -550,7 +550,7 @@ L50: goto L60; } i__2 = j + i__ * a_dim1; - if ((a[i__2].r != 0.) || (d_imag(&a[j + i__ * a_dim1]) != 0.)) { + if (a[i__2].r != 0. || d_imag(&a[j + i__ * a_dim1]) != 0.) { goto L70; } L60: @@ -581,7 +581,7 @@ L90: goto L100; } i__3 = i__ + j * a_dim1; - if ((a[i__3].r != 0.) || (d_imag(&a[i__ + j * a_dim1]) != 0.)) { + if (a[i__3].r != 0. || d_imag(&a[i__ + j * a_dim1]) != 0.) { goto L110; } L100: @@ -646,7 +646,7 @@ L150: /* Guard against zero C or R due to underflow. */ - if ((c__ == 0.) || (r__ == 0.)) { + if (c__ == 0. || r__ == 0.) { goto L200; } g = r__ / 8.; @@ -657,8 +657,7 @@ L160: d__1 = max(f,c__); /* Computing MIN */ d__2 = min(r__,g); - if (((c__ >= g) || (max(d__1,ca) >= sfmax2)) || (min(d__2,ra) <= - sfmin2)) { + if (c__ >= g || max(d__1,ca) >= sfmax2 || min(d__2,ra) <= sfmin2) { goto L170; } f *= 8.; @@ -674,8 +673,7 @@ L170: L180: /* Computing MIN */ d__1 = min(f,c__), d__1 = min(d__1,g); - if (((g < r__) || (max(r__,ra) >= sfmax2)) || (min(d__1,ca) <= sfmin2) - ) { + if (g < r__ || max(r__,ra) >= sfmax2 || min(d__1,ca) <= sfmin2) { goto L190; } f /= 8.; @@ -1569,9 +1567,9 @@ L210: *info = -3; } else if (*lda < max(1,*n)) { *info = -5; - } else if ((*ldvl < 1) || (wantvl && *ldvl < *n)) { + } else if (*ldvl < 1 || wantvl && *ldvl < *n) { *info = -8; - } else if ((*ldvr < 1) || (wantvr && *ldvr < *n)) { + } else if (*ldvr < 1 || wantvr && *ldvr < *n) { *info = -10; } @@ -1589,12 +1587,12 @@ L210: */ minwrk = 1; - if (*info == 0 && ((*lwork >= 1) || (lquery))) { + if (*info == 0 && (*lwork >= 1 || lquery)) { maxwrk = *n + *n * ilaenv_(&c__1, "ZGEHRD", " ", n, &c__1, n, &c__0, ( ftnlen)6, (ftnlen)1); if (! wantvl && ! wantvr) { /* Computing MAX */ - i__1 = 1, i__2 = (*n) << (1); + i__1 = 1, i__2 = *n << 1; minwrk = max(i__1,i__2); /* Computing MAX */ i__1 = ilaenv_(&c__8, "ZHSEQR", "EN", n, &c__1, n, &c_n1, (ftnlen) @@ -1609,12 +1607,12 @@ L210: i__1 = min(maxb,*n), i__2 = max(i__3,i__4); k = min(i__1,i__2); /* Computing MAX */ - i__1 = k * (k + 2), i__2 = (*n) << (1); + i__1 = k * (k + 2), i__2 = *n << 1; hswork = max(i__1,i__2); maxwrk = max(maxwrk,hswork); } else { /* Computing MAX */ - i__1 = 1, i__2 = (*n) << (1); + i__1 = 1, i__2 = *n << 1; minwrk = max(i__1,i__2); /* Computing MAX */ i__1 = maxwrk, i__2 = *n + (*n - 1) * ilaenv_(&c__1, "ZUNGHR", @@ -1633,10 +1631,10 @@ L210: i__1 = min(maxb,*n), i__2 = max(i__3,i__4); k = min(i__1,i__2); /* Computing MAX */ - i__1 = k * (k + 2), i__2 = (*n) << (1); + i__1 = k * (k + 2), i__2 = *n << 1; hswork = max(i__1,i__2); /* Computing MAX */ - i__1 = max(maxwrk,hswork), i__2 = (*n) << (1); + i__1 = max(maxwrk,hswork), i__2 = *n << 1; maxwrk = max(i__1,i__2); } work[1].r = (doublereal) maxwrk, work[1].i = 0.; @@ -1799,7 +1797,7 @@ L210: goto L50; } - if ((wantvl) || (wantvr)) { + if (wantvl || wantvr) { /* Compute left and/or right eigenvectors @@ -2039,9 +2037,9 @@ L50: *info = 0; if (*n < 0) { *info = -1; - } else if ((*ilo < 1) || (*ilo > max(1,*n))) { + } else if (*ilo < 1 || *ilo > max(1,*n)) { *info = -2; - } else if ((*ihi < min(*ilo,*n)) || (*ihi > *n)) { + } else if (*ihi < min(*ilo,*n) || *ihi > *n) { *info = -3; } else if (*lda < max(1,*n)) { *info = -5; @@ -2242,9 +2240,9 @@ L50: lquery = *lwork == -1; if (*n < 0) { *info = -1; - } else if ((*ilo < 1) || (*ilo > max(1,*n))) { + } else if (*ilo < 1 || *ilo > max(1,*n)) { *info = -2; - } else if ((*ihi < min(*ilo,*n)) || (*ihi > *n)) { + } else if (*ihi < min(*ilo,*n) || *ihi > *n) { *info = -3; } else if (*lda < max(1,*n)) { *info = -5; @@ -2322,7 +2320,7 @@ L50: } ldwork = *n; - if ((nb < nbmin) || (nb >= nh)) { + if (nb < nbmin || nb >= nh) { /* Use unblocked code below */ @@ -3034,24 +3032,23 @@ L50: Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + (mm + *n) * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + (mm + *n) * ilaenv_(&c__1, "ZGEBRD", " ", &mm, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1) ; maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *nrhs * ilaenv_(&c__1, - "ZUNMBR", "QLC", &mm, nrhs, n, &c_n1, (ftnlen)6, (ftnlen) - 3); + i__1 = maxwrk, i__2 = (*n << 1) + *nrhs * ilaenv_(&c__1, "ZUNMBR", + "QLC", &mm, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + (*n - 1) * ilaenv_(&c__1, - "ZUNMBR", "PLN", n, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); + i__1 = maxwrk, i__2 = (*n << 1) + (*n - 1) * ilaenv_(&c__1, "ZUN" + "MBR", "PLN", n, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * *nrhs; + i__1 = maxwrk, i__2 = (*n << 1) + *n * *nrhs; maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = ((*n) << (1)) + mm, i__2 = ((*n) << (1)) + *n * *nrhs; + i__1 = (*n << 1) + mm, i__2 = (*n << 1) + *n * *nrhs; minwrk = max(i__1,i__2); } if (*n > *m) { @@ -3065,17 +3062,17 @@ L50: maxwrk = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); /* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + ((*m) << (2)) + ((*m) << (1)) - * ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, ( + i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m << 1) * + ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, ( ftnlen)6, (ftnlen)1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + ((*m) << (2)) + *nrhs * - ilaenv_(&c__1, "ZUNMBR", "QLC", m, nrhs, m, &c_n1, ( - ftnlen)6, (ftnlen)3); + i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + *nrhs * ilaenv_(& + c__1, "ZUNMBR", "QLC", m, nrhs, m, &c_n1, (ftnlen)6, ( + ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + ((*m) << (2)) + (*m - 1) * + i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m - 1) * ilaenv_(&c__1, "ZUNMLQ", "LC", n, nrhs, m, &c_n1, ( ftnlen)6, (ftnlen)2); maxwrk = max(i__1,i__2); @@ -3085,34 +3082,33 @@ L50: maxwrk = max(i__1,i__2); } else { /* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + ((*m) << (1)); + i__1 = maxwrk, i__2 = *m * *m + (*m << 1); maxwrk = max(i__1,i__2); } /* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + ((*m) << (2)) + *m * *nrhs; + i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + *m * *nrhs; maxwrk = max(i__1,i__2); } else { /* Path 2 - underdetermined. */ - maxwrk = ((*m) << (1)) + (*n + *m) * ilaenv_(&c__1, "ZGEBRD", + maxwrk = (*m << 1) + (*n + *m) * ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *nrhs * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *nrhs * ilaenv_(&c__1, "ZUNMBR", "QLC", m, nrhs, m, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, - "ZUNMBR", "PLN", n, nrhs, m, &c_n1, (ftnlen)6, ( - ftnlen)3); + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR" + , "PLN", n, nrhs, m, &c_n1, (ftnlen)6, (ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * *nrhs; + i__1 = maxwrk, i__2 = (*m << 1) + *m * *nrhs; maxwrk = max(i__1,i__2); } /* Computing MAX */ - i__1 = ((*m) << (1)) + *n, i__2 = ((*m) << (1)) + *m * *nrhs; + i__1 = (*m << 1) + *n, i__2 = (*m << 1) + *m * *nrhs; minwrk = max(i__1,i__2); } minwrk = min(minwrk,maxwrk); @@ -3134,7 +3130,7 @@ L50: /* Quick return if possible. */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { *rank = 0; return 0; } @@ -3290,10 +3286,9 @@ L50: } else /* if(complicated condition) */ { /* Computing MAX */ - i__1 = *m, i__2 = ((*m) << (1)) - 4, i__1 = max(i__1,i__2), i__1 = - max(i__1,*nrhs), i__2 = *n - *m * 3; - if (*n >= mnthr && *lwork >= ((*m) << (2)) + *m * *m + max(i__1,i__2)) - { + i__1 = *m, i__2 = (*m << 1) - 4, i__1 = max(i__1,i__2), i__1 = max( + i__1,*nrhs), i__2 = *n - *m * 3; + if (*n >= mnthr && *lwork >= (*m << 2) + *m * *m + max(i__1,i__2)) { /* Path 2a - underdetermined, with many more columns than rows @@ -3305,10 +3300,10 @@ L50: Computing MAX Computing MAX */ - i__3 = *m, i__4 = ((*m) << (1)) - 4, i__3 = max(i__3,i__4), i__3 = - max(i__3,*nrhs), i__4 = *n - *m * 3; - i__1 = ((*m) << (2)) + *m * *lda + max(i__3,i__4), i__2 = *m * * - lda + *m + *m * *nrhs; + i__3 = *m, i__4 = (*m << 1) - 4, i__3 = max(i__3,i__4), i__3 = + max(i__3,*nrhs), i__4 = *n - *m * 3; + i__1 = (*m << 2) + *m * *lda + max(i__3,i__4), i__2 = *m * *lda + + *m + *m * *nrhs; if (*lwork >= max(i__1,i__2)) { ldwork = *lda; } @@ -4084,14 +4079,14 @@ L10: mnthr2 = (integer) (minmn * 5. / 3.); wntqa = lsame_(jobz, "A"); wntqs = lsame_(jobz, "S"); - wntqas = (wntqa) || (wntqs); + wntqas = wntqa || wntqs; wntqo = lsame_(jobz, "O"); wntqn = lsame_(jobz, "N"); minwrk = 1; maxwrk = 1; lquery = *lwork == -1; - if (! ((((wntqa) || (wntqs)) || (wntqo)) || (wntqn))) { + if (! (wntqa || wntqs || wntqo || wntqn)) { *info = -1; } else if (*m < 0) { *info = -2; @@ -4099,11 +4094,11 @@ L10: *info = -3; } else if (*lda < max(1,*m)) { *info = -5; - } else if (((*ldu < 1) || (wntqas && *ldu < *m)) || (wntqo && *m < *n && * - ldu < *m)) { + } else if (*ldu < 1 || wntqas && *ldu < *m || wntqo && *m < *n && *ldu < * + m) { *info = -8; - } else if ((((*ldvt < 1) || (wntqa && *ldvt < *n)) || (wntqs && *ldvt < - minmn)) || (wntqo && *m >= *n && *ldvt < *n)) { + } else if (*ldvt < 1 || wntqa && *ldvt < *n || wntqs && *ldvt < minmn || + wntqo && *m >= *n && *ldvt < *n) { *info = -10; } @@ -4134,9 +4129,9 @@ L10: wrkbl = *n + *n * ilaenv_(&c__1, "ZGEQRF", " ", m, n, & c_n1, &c_n1, (ftnlen)6, (ftnlen)1); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + ((*n) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); maxwrk = wrkbl; minwrk = *n * 3; @@ -4151,22 +4146,22 @@ L10: " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + ((*n) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); maxwrk = *m * *n + *n * *n + wrkbl; - minwrk = ((*n) << (1)) * *n + *n * 3; + minwrk = (*n << 1) * *n + *n * 3; } else if (wntqs) { /* Path 3 (M much larger than N, JOBZ='S') */ @@ -4178,17 +4173,17 @@ L10: " ", m, n, n, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + ((*n) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); @@ -4205,38 +4200,38 @@ L10: " ", m, m, n, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + ((*n) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*n << 1) + (*n << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", n, n, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "QLN", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); maxwrk = *n * *n + wrkbl; - minwrk = *n * *n + ((*n) << (1)) + *m; + minwrk = *n * *n + (*n << 1) + *m; } } else if (*m >= mnthr2) { /* Path 5 (M much larger than N, but not as much as MNTHR1) */ - maxwrk = ((*n) << (1)) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", + maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - minwrk = ((*n) << (1)) + *m; + minwrk = (*n << 1) + *m; if (wntqo) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "Q", m, n, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); @@ -4244,23 +4239,23 @@ L10: minwrk += *n * *n; } else if (wntqs) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "Q", m, n, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); } else if (wntqa) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "P", n, n, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); @@ -4269,17 +4264,17 @@ L10: /* Path 6 (M at least N, but not much larger) */ - maxwrk = ((*n) << (1)) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", + maxwrk = (*n << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - minwrk = ((*n) << (1)) + *m; + minwrk = (*n << 1) + *m; if (wntqo) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "QLN", m, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); @@ -4287,23 +4282,23 @@ L10: minwrk += *n * *n; } else if (wntqs) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNMBR", "QLN", m, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); } else if (wntqa) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "PRC", n, n, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*n) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*n << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "QLN", m, m, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); @@ -4325,9 +4320,9 @@ L10: maxwrk = *m + *m * ilaenv_(&c__1, "ZGELQF", " ", m, n, & c_n1, &c_n1, (ftnlen)6, (ftnlen)1); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + ((*m) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = maxwrk, i__2 = (*m << 1) + (*m << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); maxwrk = max(i__1,i__2); minwrk = *m * 3; } else if (wntqo) { @@ -4341,22 +4336,22 @@ L10: " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + ((*m) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "PRC", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "QLN", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); maxwrk = *m * *n + *m * *m + wrkbl; - minwrk = ((*m) << (1)) * *m + *m * 3; + minwrk = (*m << 1) * *m + *m * 3; } else if (wntqs) { /* Path 3t (N much larger than M, JOBZ='S') */ @@ -4368,17 +4363,17 @@ L10: " ", m, n, m, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + ((*m) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "PRC", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "QLN", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); @@ -4395,38 +4390,38 @@ L10: " ", n, n, m, &c_n1, (ftnlen)6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + ((*m) << (1)) * - ilaenv_(&c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, - (ftnlen)6, (ftnlen)1); + i__1 = wrkbl, i__2 = (*m << 1) + (*m << 1) * ilaenv_(& + c__1, "ZGEBRD", " ", m, m, &c_n1, &c_n1, (ftnlen) + 6, (ftnlen)1); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "PRC", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); /* Computing MAX */ - i__1 = wrkbl, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = wrkbl, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "QLN", m, m, m, &c_n1, (ftnlen)6, ( ftnlen)3); wrkbl = max(i__1,i__2); maxwrk = *m * *m + wrkbl; - minwrk = *m * *m + ((*m) << (1)) + *n; + minwrk = *m * *m + (*m << 1) + *n; } } else if (*n >= mnthr2) { /* Path 5t (N much larger than M, but not as much as MNTHR1) */ - maxwrk = ((*m) << (1)) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", + maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - minwrk = ((*m) << (1)) + *n; + minwrk = (*m << 1) + *n; if (wntqo) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "P", m, n, m, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); @@ -4434,23 +4429,23 @@ L10: minwrk += *m * *m; } else if (wntqs) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "P", m, n, m, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); } else if (wntqa) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "P", n, n, m, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "Q", m, m, n, &c_n1, (ftnlen)6, (ftnlen) 1); maxwrk = max(i__1,i__2); @@ -4459,17 +4454,17 @@ L10: /* Path 6t (N greater than M, but not much larger) */ - maxwrk = ((*m) << (1)) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", + maxwrk = (*m << 1) + (*m + *n) * ilaenv_(&c__1, "ZGEBRD", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - minwrk = ((*m) << (1)) + *n; + minwrk = (*m << 1) + *n; if (wntqo) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "PRC", m, n, m, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNMBR", "QLN", m, m, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); @@ -4477,23 +4472,23 @@ L10: minwrk += *m * *m; } else if (wntqs) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "PRC", m, n, m, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "QLN", m, m, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); } else if (wntqa) { /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *n * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *n * ilaenv_(&c__1, "ZUNGBR", "PRC", n, n, m, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); /* Computing MAX */ - i__1 = maxwrk, i__2 = ((*m) << (1)) + *m * ilaenv_(&c__1, + i__1 = maxwrk, i__2 = (*m << 1) + *m * ilaenv_(&c__1, "ZUNGBR", "QLN", m, m, n, &c_n1, (ftnlen)6, ( ftnlen)3); maxwrk = max(i__1,i__2); @@ -4517,7 +4512,7 @@ L10: /* Quick return if possible */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { if (*lwork >= 1) { work[1].r = 1., work[1].i = 0.; } @@ -5462,8 +5457,8 @@ L10: i__2 = *m - 1; i__1 = *m - 1; - zlaset_("U", &i__2, &i__1, &c_b59, &c_b59, &a[((a_dim1) << (1) - ) + 1], lda); + zlaset_("U", &i__2, &i__1, &c_b59, &c_b59, &a[(a_dim1 << 1) + + 1], lda); ie = 1; itauq = 1; itaup = itauq + *m; @@ -5773,8 +5768,8 @@ L10: i__1 = *m - 1; i__2 = *m - 1; - zlaset_("U", &i__1, &i__2, &c_b59, &c_b59, &a[((a_dim1) << (1) - ) + 1], lda); + zlaset_("U", &i__1, &i__2, &c_b59, &c_b59, &a[(a_dim1 << 1) + + 1], lda); ie = 1; itauq = itau; itaup = itauq + *m; @@ -6562,7 +6557,7 @@ L10: /* Quick return if possible */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } @@ -6575,7 +6570,7 @@ L10: jp = j - 1 + izamax_(&i__2, &a[j + j * a_dim1], &c__1); ipiv[j] = jp; i__2 = jp + j * a_dim1; - if ((a[i__2].r != 0.) || (a[i__2].i != 0.)) { + if (a[i__2].r != 0. || a[i__2].i != 0.) { /* Apply the interchange to columns 1:N. */ @@ -6717,7 +6712,7 @@ L10: /* Quick return if possible */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } @@ -6725,7 +6720,7 @@ L10: nb = ilaenv_(&c__1, "ZGETRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen) 1); - if ((nb <= 1) || (nb >= min(*m,*n))) { + if (nb <= 1 || nb >= min(*m,*n)) { /* Use unblocked code. */ @@ -6914,7 +6909,7 @@ L10: /* Quick return if possible */ - if ((*n == 0) || (*nrhs == 0)) { + if (*n == 0 || *nrhs == 0) { return 0; } @@ -7150,7 +7145,7 @@ L10: /* Function Body */ wantz = lsame_(jobz, "V"); lower = lsame_(uplo, "L"); - lquery = ((*lwork == -1) || (*lrwork == -1)) || (*liwork == -1); + lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1; *info = 0; if (*n <= 1) { @@ -7162,10 +7157,10 @@ L10: liopt = liwmin; } else { if (wantz) { - lwmin = ((*n) << (1)) + *n * *n; + lwmin = (*n << 1) + *n * *n; /* Computing 2nd power */ i__1 = *n; - lrwmin = *n * 5 + 1 + ((i__1 * i__1) << (1)); + lrwmin = *n * 5 + 1 + (i__1 * i__1 << 1); liwmin = *n * 5 + 3; } else { lwmin = *n + 1; @@ -7176,9 +7171,9 @@ L10: lropt = lrwmin; liopt = liwmin; } - if (! ((wantz) || (lsame_(jobz, "N")))) { + if (! (wantz || lsame_(jobz, "N"))) { *info = -1; - } else if (! ((lower) || (lsame_(uplo, "U")))) { + } else if (! (lower || lsame_(uplo, "U"))) { *info = -2; } else if (*n < 0) { *info = -3; @@ -7505,7 +7500,7 @@ L10: i__1 = i__; e[i__1] = alpha.r; - if ((taui.r != 0.) || (taui.i != 0.)) { + if (taui.r != 0. || taui.i != 0.) { /* Apply H(i) from both sides to A(1:i,1:i) */ @@ -7584,7 +7579,7 @@ L10: i__2 = i__; e[i__2] = alpha.r; - if ((taui.r != 0.) || (taui.i != 0.)) { + if (taui.r != 0. || taui.i != 0.) { /* Apply H(i) from both sides to A(i+1:n,i+1:n) */ @@ -8188,7 +8183,7 @@ L10: /* Function Body */ wantt = lsame_(job, "S"); initz = lsame_(compz, "I"); - wantz = (initz) || (lsame_(compz, "V")); + wantz = initz || lsame_(compz, "V"); *info = 0; i__1 = max(1,*n); @@ -8200,13 +8195,13 @@ L10: *info = -2; } else if (*n < 0) { *info = -3; - } else if ((*ilo < 1) || (*ilo > max(1,*n))) { + } else if (*ilo < 1 || *ilo > max(1,*n)) { *info = -4; - } else if ((*ihi < min(*ilo,*n)) || (*ihi > *n)) { + } else if (*ihi < min(*ilo,*n) || *ihi > *n) { *info = -5; } else if (*ldh < max(1,*n)) { *info = -7; - } else if ((*ldz < 1) || (wantz && *ldz < max(1,*n))) { + } else if (*ldz < 1 || wantz && *ldz < max(1,*n)) { *info = -10; } else if (*lwork < max(1,*n) && ! lquery) { *info = -12; @@ -8336,7 +8331,7 @@ L10: s_cat(ch__1, a__1, i__4, &c__2, (ftnlen)2); maxb = ilaenv_(&c__8, "ZHSEQR", ch__1, n, ilo, ihi, &c_n1, (ftnlen)6, ( ftnlen)2); - if (((ns <= 1) || (ns > nh)) || (maxb >= nh)) { + if (ns <= 1 || ns > nh || maxb >= nh) { /* Use the standard double-shift algorithm */ @@ -8439,7 +8434,7 @@ L80: i2 = i__; } - if ((its == 20) || (its == 30)) { + if (its == 20 || its == 30) { /* Exceptional shifts. */ @@ -8828,7 +8823,7 @@ L180: y -= y_offset; /* Function Body */ - if ((*m <= 0) || (*n <= 0)) { + if (*m <= 0 || *n <= 0) { return 0; } @@ -9516,7 +9511,7 @@ L180: --rwork; /* Function Body */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } @@ -9806,7 +9801,7 @@ L10: if (iwork[subpbs] > smlsiz) { for (j = subpbs; j >= 1; --j) { iwork[j * 2] = (iwork[j] + 1) / 2; - iwork[((j) << (1)) - 1] = iwork[j] / 2; + iwork[(j << 1) - 1] = iwork[j] / 2; /* L20: */ } ++tlvls; @@ -9834,7 +9829,7 @@ L10: /* L40: */ } - indxq = ((*n) << (2)) + 3; + indxq = (*n << 2) + 3; /* Set up workspaces for eigenvalues only/accumulate new vectors @@ -9856,7 +9851,7 @@ L10: igivcl = igivpt + *n * lgn; igivnm = 1; - iq = igivnm + ((*n) << (1)) * lgn; + iq = igivnm + (*n << 1) * lgn; /* Computing 2nd power */ i__1 = *n; iwrem = iq + i__1 * i__1 + 1; @@ -10044,7 +10039,7 @@ L80: The first stage consists of deflating the size of the problem when there are multiple eigenvalues or if there is a zero in - the Z vector. For each such occurrence the dimension of the + the Z vector. For each such occurence the dimension of the secular equation problem is reduced by one. This stage is performed by the routine DLAED2. @@ -10184,7 +10179,7 @@ L80: */ if (*n < 0) { *info = -1; - } else if ((min(1,*n) > *cutpnt) || (*n < *cutpnt)) { + } else if (min(1,*n) > *cutpnt || *n < *cutpnt) { *info = -2; } else if (*qsiz < *n) { *info = -3; @@ -10253,8 +10248,7 @@ L80: zlaed8_(&k, n, qsiz, &q[q_offset], ldq, &d__[1], rho, cutpnt, &rwork[iz], &rwork[idlmda], &work[1], qsiz, &rwork[iw], &iwork[indxp], &iwork[ indx], &indxq[1], &perm[prmptr[curr]], &givptr[curr + 1], &givcol[ - ((givptr[curr]) << (1)) + 1], &givnum[((givptr[curr]) << (1)) + 1] - , info); + (givptr[curr] << 1) + 1], &givnum[(givptr[curr] << 1) + 1], info); prmptr[curr + 1] = prmptr[curr] + *n; givptr[curr + 1] += givptr[curr]; @@ -10478,7 +10472,7 @@ L80: *info = -3; } else if (*ldq < max(1,*n)) { *info = -5; - } else if ((*cutpnt < min(1,*n)) || (*cutpnt > *n)) { + } else if (*cutpnt < min(1,*n) || *cutpnt > *n) { *info = -8; } else if (*ldq2 < max(1,*n)) { *info = -12; @@ -10628,10 +10622,10 @@ L70: /* Record the appropriate Givens rotation */ ++(*givptr); - givcol[((*givptr) << (1)) + 1] = indxq[indx[jlam]]; - givcol[((*givptr) << (1)) + 2] = indxq[indx[j]]; - givnum[((*givptr) << (1)) + 1] = c__; - givnum[((*givptr) << (1)) + 2] = s; + givcol[(*givptr << 1) + 1] = indxq[indx[jlam]]; + givcol[(*givptr << 1) + 2] = indxq[indx[j]]; + givnum[(*givptr << 1) + 1] = c__; + givnum[(*givptr << 1) + 2] = s; zdrot_(qsiz, &q[indxq[indx[jlam]] * q_dim1 + 1], &c__1, &q[indxq[ indx[j]] * q_dim1 + 1], &c__1, &c__, &s); t = d__[jlam] * c__ * c__ + d__[j] * s * s; @@ -10962,7 +10956,7 @@ L30: i2 = i__; } - if ((its == 10) || (its == 20)) { + if (its == 10 || its == 20) { /* Exceptional shift. */ @@ -10982,7 +10976,7 @@ L30: d__1 = h__[i__3].r; z__1.r = d__1 * h__[i__2].r, z__1.i = d__1 * h__[i__2].i; u.r = z__1.r, u.i = z__1.i; - if ((u.r != 0.) || (u.i != 0.)) { + if (u.r != 0. || u.i != 0.) { i__2 = i__ - 1 + (i__ - 1) * h_dim1; z__2.r = h__[i__2].r - t.r, z__2.i = h__[i__2].i - t.i; z__1.r = z__2.r * .5, z__1.i = z__2.i * .5; @@ -11784,13 +11778,13 @@ L130: /* Function Body */ *info = 0; - if ((*icompq < 0) || (*icompq > 1)) { + if (*icompq < 0 || *icompq > 1) { *info = -1; } else if (*nl < 1) { *info = -2; } else if (*nr < 1) { *info = -3; - } else if ((*sqre < 0) || (*sqre > 1)) { + } else if (*sqre < 0 || *sqre > 1) { *info = -4; } @@ -11830,10 +11824,9 @@ L130: i__1 = *givptr; for (i__ = 1; i__ <= i__1; ++i__) { - zdrot_(nrhs, &b[givcol[i__ + ((givcol_dim1) << (1))] + b_dim1], - ldb, &b[givcol[i__ + givcol_dim1] + b_dim1], ldb, &givnum[ - i__ + ((givnum_dim1) << (1))], &givnum[i__ + givnum_dim1]) - ; + zdrot_(nrhs, &b[givcol[i__ + (givcol_dim1 << 1)] + b_dim1], ldb, & + b[givcol[i__ + givcol_dim1] + b_dim1], ldb, &givnum[i__ + + (givnum_dim1 << 1)], &givnum[i__ + givnum_dim1]); /* L10: */ } @@ -11862,41 +11855,40 @@ L130: for (j = 1; j <= i__1; ++j) { diflj = difl[j]; dj = poles[j + poles_dim1]; - dsigj = -poles[j + ((poles_dim1) << (1))]; + dsigj = -poles[j + (poles_dim1 << 1)]; if (j < *k) { difrj = -difr[j + difr_dim1]; - dsigjp = -poles[j + 1 + ((poles_dim1) << (1))]; + dsigjp = -poles[j + 1 + (poles_dim1 << 1)]; } - if ((z__[j] == 0.) || (poles[j + ((poles_dim1) << (1))] == 0.) - ) { + if (z__[j] == 0. || poles[j + (poles_dim1 << 1)] == 0.) { rwork[j] = 0.; } else { - rwork[j] = -poles[j + ((poles_dim1) << (1))] * z__[j] / - diflj / (poles[j + ((poles_dim1) << (1))] + dj); + rwork[j] = -poles[j + (poles_dim1 << 1)] * z__[j] / diflj + / (poles[j + (poles_dim1 << 1)] + dj); } i__2 = j - 1; for (i__ = 1; i__ <= i__2; ++i__) { - if ((z__[i__] == 0.) || (poles[i__ + ((poles_dim1) << (1)) - ] == 0.)) { + if (z__[i__] == 0. || poles[i__ + (poles_dim1 << 1)] == + 0.) { rwork[i__] = 0.; } else { - rwork[i__] = poles[i__ + ((poles_dim1) << (1))] * z__[ - i__] / (dlamc3_(&poles[i__ + ((poles_dim1) << - (1))], &dsigj) - diflj) / (poles[i__ + (( - poles_dim1) << (1))] + dj); + rwork[i__] = poles[i__ + (poles_dim1 << 1)] * z__[i__] + / (dlamc3_(&poles[i__ + (poles_dim1 << 1)], & + dsigj) - diflj) / (poles[i__ + (poles_dim1 << + 1)] + dj); } /* L30: */ } i__2 = *k; for (i__ = j + 1; i__ <= i__2; ++i__) { - if ((z__[i__] == 0.) || (poles[i__ + ((poles_dim1) << (1)) - ] == 0.)) { + if (z__[i__] == 0. || poles[i__ + (poles_dim1 << 1)] == + 0.) { rwork[i__] = 0.; } else { - rwork[i__] = poles[i__ + ((poles_dim1) << (1))] * z__[ - i__] / (dlamc3_(&poles[i__ + ((poles_dim1) << - (1))], &dsigjp) + difrj) / (poles[i__ + (( - poles_dim1) << (1))] + dj); + rwork[i__] = poles[i__ + (poles_dim1 << 1)] * z__[i__] + / (dlamc3_(&poles[i__ + (poles_dim1 << 1)], & + dsigjp) + difrj) / (poles[i__ + (poles_dim1 << + 1)] + dj); } /* L40: */ } @@ -11911,7 +11903,7 @@ L130: $ B( J, 1 ), LDB ) */ - i__ = *k + ((*nrhs) << (1)); + i__ = *k + (*nrhs << 1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = *k; @@ -11923,10 +11915,9 @@ L130: } /* L60: */ } - dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + ((*nrhs) << (1) - )], k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1], & - c__1); - i__ = *k + ((*nrhs) << (1)); + dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + (*nrhs << 1)], + k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1], &c__1); + i__ = *k + (*nrhs << 1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = *k; @@ -11937,9 +11928,9 @@ L130: } /* L80: */ } - dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + ((*nrhs) << (1) - )], k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1 + * - nrhs], &c__1); + dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + (*nrhs << 1)], + k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1 + *nrhs], + &c__1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = j + jcol * b_dim1; @@ -11976,23 +11967,22 @@ L130: } else { i__1 = *k; for (j = 1; j <= i__1; ++j) { - dsigj = poles[j + ((poles_dim1) << (1))]; + dsigj = poles[j + (poles_dim1 << 1)]; if (z__[j] == 0.) { rwork[j] = 0.; } else { rwork[j] = -z__[j] / difl[j] / (dsigj + poles[j + - poles_dim1]) / difr[j + ((difr_dim1) << (1))]; + poles_dim1]) / difr[j + (difr_dim1 << 1)]; } i__2 = j - 1; for (i__ = 1; i__ <= i__2; ++i__) { if (z__[j] == 0.) { rwork[i__] = 0.; } else { - d__1 = -poles[i__ + 1 + ((poles_dim1) << (1))]; + d__1 = -poles[i__ + 1 + (poles_dim1 << 1)]; rwork[i__] = z__[j] / (dlamc3_(&dsigj, &d__1) - difr[ i__ + difr_dim1]) / (dsigj + poles[i__ + - poles_dim1]) / difr[i__ + ((difr_dim1) << (1)) - ]; + poles_dim1]) / difr[i__ + (difr_dim1 << 1)]; } /* L110: */ } @@ -12001,10 +11991,10 @@ L130: if (z__[j] == 0.) { rwork[i__] = 0.; } else { - d__1 = -poles[i__ + ((poles_dim1) << (1))]; + d__1 = -poles[i__ + (poles_dim1 << 1)]; rwork[i__] = z__[j] / (dlamc3_(&dsigj, &d__1) - difl[ i__]) / (dsigj + poles[i__ + poles_dim1]) / - difr[i__ + ((difr_dim1) << (1))]; + difr[i__ + (difr_dim1 << 1)]; } /* L120: */ } @@ -12017,7 +12007,7 @@ L130: $ BX( J, 1 ), LDBX ) */ - i__ = *k + ((*nrhs) << (1)); + i__ = *k + (*nrhs << 1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = *k; @@ -12029,10 +12019,9 @@ L130: } /* L140: */ } - dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + ((*nrhs) << (1) - )], k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1], & - c__1); - i__ = *k + ((*nrhs) << (1)); + dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + (*nrhs << 1)], + k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1], &c__1); + i__ = *k + (*nrhs << 1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = *k; @@ -12043,9 +12032,9 @@ L130: } /* L160: */ } - dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + ((*nrhs) << (1) - )], k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1 + * - nrhs], &c__1); + dgemv_("T", k, nrhs, &c_b1015, &rwork[*k + 1 + (*nrhs << 1)], + k, &rwork[1], &c__1, &c_b324, &rwork[*k + 1 + *nrhs], + &c__1); i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = j + jcol * bx_dim1; @@ -12092,9 +12081,9 @@ L130: for (i__ = *givptr; i__ >= 1; --i__) { d__1 = -givnum[i__ + givnum_dim1]; - zdrot_(nrhs, &b[givcol[i__ + ((givcol_dim1) << (1))] + b_dim1], - ldb, &b[givcol[i__ + givcol_dim1] + b_dim1], ldb, &givnum[ - i__ + ((givnum_dim1) << (1))], &d__1); + zdrot_(nrhs, &b[givcol[i__ + (givcol_dim1 << 1)] + b_dim1], ldb, & + b[givcol[i__ + givcol_dim1] + b_dim1], ldb, &givnum[i__ + + (givnum_dim1 << 1)], &d__1); /* L200: */ } } @@ -12333,7 +12322,7 @@ L130: /* Function Body */ *info = 0; - if ((*icompq < 0) || (*icompq > 1)) { + if (*icompq < 0 || *icompq > 1) { *info = -1; } else if (*smlsiz < 3) { *info = -2; @@ -12408,7 +12397,7 @@ L130: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) */ - j = (nl * *nrhs) << (1); + j = nl * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nlf + nl - 1; @@ -12421,9 +12410,8 @@ L130: /* L20: */ } dgemm_("T", "N", &nl, nrhs, &nl, &c_b1015, &u[nlf + u_dim1], ldu, & - rwork[((nl * *nrhs) << (1)) + 1], &nl, &c_b324, &rwork[1], & - nl); - j = (nl * *nrhs) << (1); + rwork[(nl * *nrhs << 1) + 1], &nl, &c_b324, &rwork[1], &nl); + j = nl * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nlf + nl - 1; @@ -12435,8 +12423,8 @@ L130: /* L40: */ } dgemm_("T", "N", &nl, nrhs, &nl, &c_b1015, &u[nlf + u_dim1], ldu, & - rwork[((nl * *nrhs) << (1)) + 1], &nl, &c_b324, &rwork[nl * * - nrhs + 1], &nl); + rwork[(nl * *nrhs << 1) + 1], &nl, &c_b324, &rwork[nl * *nrhs + + 1], &nl); jreal = 0; jimag = nl * *nrhs; i__2 = *nrhs; @@ -12463,7 +12451,7 @@ L130: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) */ - j = (nr * *nrhs) << (1); + j = nr * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nrf + nr - 1; @@ -12476,9 +12464,8 @@ L130: /* L80: */ } dgemm_("T", "N", &nr, nrhs, &nr, &c_b1015, &u[nrf + u_dim1], ldu, & - rwork[((nr * *nrhs) << (1)) + 1], &nr, &c_b324, &rwork[1], & - nr); - j = (nr * *nrhs) << (1); + rwork[(nr * *nrhs << 1) + 1], &nr, &c_b324, &rwork[1], &nr); + j = nr * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nrf + nr - 1; @@ -12490,8 +12477,8 @@ L130: /* L100: */ } dgemm_("T", "N", &nr, nrhs, &nr, &c_b1015, &u[nrf + u_dim1], ldu, & - rwork[((nr * *nrhs) << (1)) + 1], &nr, &c_b324, &rwork[nr * * - nrhs + 1], &nr); + rwork[(nr * *nrhs << 1) + 1], &nr, &c_b324, &rwork[nr * *nrhs + + 1], &nr); jreal = 0; jimag = nr * *nrhs; i__2 = *nrhs; @@ -12534,7 +12521,7 @@ L130: sqre = 0; for (lvl = nlvl; lvl >= 1; --lvl) { - lvl2 = ((lvl) << (1)) - 1; + lvl2 = (lvl << 1) - 1; /* find the first node LF and last node LL on @@ -12547,7 +12534,7 @@ L130: } else { i__1 = lvl - 1; lf = pow_ii(&c__2, &i__1); - ll = ((lf) << (1)) - 1; + ll = (lf << 1) - 1; } i__1 = ll; for (i__ = lf; i__ <= i__1; ++i__) { @@ -12583,7 +12570,7 @@ L170: j = 0; i__1 = nlvl; for (lvl = 1; lvl <= i__1; ++lvl) { - lvl2 = ((lvl) << (1)) - 1; + lvl2 = (lvl << 1) - 1; /* Find the first node LF and last node LL on @@ -12596,7 +12583,7 @@ L170: } else { i__2 = lvl - 1; lf = pow_ii(&c__2, &i__2); - ll = ((lf) << (1)) - 1; + ll = (lf << 1) - 1; } i__2 = lf; for (i__ = ll; i__ >= i__2; --i__) { @@ -12654,7 +12641,7 @@ L170: $ B( NLF, 1 ), LDB, ZERO, BX( NLF, 1 ), LDBX ) */ - j = (nlp1 * *nrhs) << (1); + j = nlp1 * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nlf + nlp1 - 1; @@ -12667,9 +12654,9 @@ L170: /* L210: */ } dgemm_("T", "N", &nlp1, nrhs, &nlp1, &c_b1015, &vt[nlf + vt_dim1], - ldu, &rwork[((nlp1 * *nrhs) << (1)) + 1], &nlp1, &c_b324, & - rwork[1], &nlp1); - j = (nlp1 * *nrhs) << (1); + ldu, &rwork[(nlp1 * *nrhs << 1) + 1], &nlp1, &c_b324, &rwork[ + 1], &nlp1); + j = nlp1 * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nlf + nlp1 - 1; @@ -12681,8 +12668,8 @@ L170: /* L230: */ } dgemm_("T", "N", &nlp1, nrhs, &nlp1, &c_b1015, &vt[nlf + vt_dim1], - ldu, &rwork[((nlp1 * *nrhs) << (1)) + 1], &nlp1, &c_b324, & - rwork[nlp1 * *nrhs + 1], &nlp1); + ldu, &rwork[(nlp1 * *nrhs << 1) + 1], &nlp1, &c_b324, &rwork[ + nlp1 * *nrhs + 1], &nlp1); jreal = 0; jimag = nlp1 * *nrhs; i__2 = *nrhs; @@ -12709,7 +12696,7 @@ L170: $ B( NRF, 1 ), LDB, ZERO, BX( NRF, 1 ), LDBX ) */ - j = (nrp1 * *nrhs) << (1); + j = nrp1 * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nrf + nrp1 - 1; @@ -12722,9 +12709,9 @@ L170: /* L270: */ } dgemm_("T", "N", &nrp1, nrhs, &nrp1, &c_b1015, &vt[nrf + vt_dim1], - ldu, &rwork[((nrp1 * *nrhs) << (1)) + 1], &nrp1, &c_b324, & - rwork[1], &nrp1); - j = (nrp1 * *nrhs) << (1); + ldu, &rwork[(nrp1 * *nrhs << 1) + 1], &nrp1, &c_b324, &rwork[ + 1], &nrp1); + j = nrp1 * *nrhs << 1; i__2 = *nrhs; for (jcol = 1; jcol <= i__2; ++jcol) { i__3 = nrf + nrp1 - 1; @@ -12736,8 +12723,8 @@ L170: /* L290: */ } dgemm_("T", "N", &nrp1, nrhs, &nrp1, &c_b1015, &vt[nrf + vt_dim1], - ldu, &rwork[((nrp1 * *nrhs) << (1)) + 1], &nrp1, &c_b324, & - rwork[nrp1 * *nrhs + 1], &nrp1); + ldu, &rwork[(nrp1 * *nrhs << 1) + 1], &nrp1, &c_b324, &rwork[ + nrp1 * *nrhs + 1], &nrp1); jreal = 0; jimag = nrp1 * *nrhs; i__2 = *nrhs; @@ -12961,7 +12948,7 @@ L330: *info = -3; } else if (*nrhs < 1) { *info = -4; - } else if ((*ldb < 1) || (*ldb < *n)) { + } else if (*ldb < 1 || *ldb < *n) { *info = -8; } if (*info != 0) { @@ -12974,7 +12961,7 @@ L330: /* Set up the tolerance. */ - if ((*rcond <= 0.) || (*rcond >= 1.)) { + if (*rcond <= 0. || *rcond >= 1.) { *rcond = eps; } @@ -13009,7 +12996,7 @@ L330: zdrot_(&c__1, &b[i__ + b_dim1], &c__1, &b[i__ + 1 + b_dim1], & c__1, &cs, &sn); } else { - rwork[((i__) << (1)) - 1] = cs; + rwork[(i__ << 1) - 1] = cs; rwork[i__ * 2] = sn; } /* L10: */ @@ -13019,7 +13006,7 @@ L330: for (i__ = 1; i__ <= i__1; ++i__) { i__2 = *n - 1; for (j = 1; j <= i__2; ++j) { - cs = rwork[((j) << (1)) - 1]; + cs = rwork[(j << 1) - 1]; sn = rwork[j * 2]; zdrot_(&c__1, &b[j + i__ * b_dim1], &c__1, &b[j + 1 + i__ * b_dim1], &c__1, &cs, &sn); @@ -13204,12 +13191,12 @@ L330: vt = *smlsiz * *n + 1; difl = vt + smlszp * *n; difr = difl + nlvl * *n; - z__ = difr + ((nlvl * *n) << (1)); + z__ = difr + (nlvl * *n << 1); c__ = z__ + nlvl * *n; s = c__ + *n; poles = s + *n; - givnum = poles + ((nlvl) << (1)) * *n; - nrwork = givnum + ((nlvl) << (1)) * *n; + givnum = poles + (nlvl << 1) * *n; + nrwork = givnum + (nlvl << 1) * *n; bx = 1; irwrb = nrwork; @@ -13221,7 +13208,7 @@ L330: givptr = k + *n; perm = givptr + *n; givcol = perm + nlvl * *n; - iwk = givcol + ((nlvl * *n) << (1)); + iwk = givcol + (nlvl * *n << 1); st = 1; sqre = 0; @@ -13239,7 +13226,7 @@ L330: i__1 = nm1; for (i__ = 1; i__ <= i__1; ++i__) { - if (((d__1 = e[i__], abs(d__1)) < eps) || (i__ == nm1)) { + if ((d__1 = e[i__], abs(d__1)) < eps || i__ == nm1) { ++nsub; iwork[nsub] = st; @@ -13612,8 +13599,8 @@ doublereal zlange_(char *norm, integer *m, integer *n, doublecomplex *a, } /* L20: */ } - } else if ((lsame_(norm, "O")) || (*(unsigned char * - )norm == '1')) { + } else if (lsame_(norm, "O") || *(unsigned char *) + norm == '1') { /* Find norm1(A). */ @@ -13655,8 +13642,7 @@ doublereal zlange_(char *norm, integer *m, integer *n, doublecomplex *a, value = max(d__1,d__2); /* L80: */ } - } else if ((lsame_(norm, "F")) || (lsame_(norm, - "E"))) { + } else if (lsame_(norm, "F") || lsame_(norm, "E")) { /* Find normF(A). */ @@ -13813,8 +13799,7 @@ doublereal zlanhe_(char *norm, char *uplo, integer *n, doublecomplex *a, /* L40: */ } } - } else if (((lsame_(norm, "I")) || (lsame_(norm, - "O"))) || (*(unsigned char *)norm == '1')) { + } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') { /* Find normI(A) ( = norm1(A), since A is hermitian). */ @@ -13862,8 +13847,7 @@ doublereal zlanhe_(char *norm, char *uplo, integer *n, doublecomplex *a, /* L100: */ } } - } else if ((lsame_(norm, "F")) || (lsame_(norm, - "E"))) { + } else if (lsame_(norm, "F") || lsame_(norm, "E")) { /* Find normF(A). */ @@ -14018,8 +14002,8 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, } /* L20: */ } - } else if ((lsame_(norm, "O")) || (*(unsigned char * - )norm == '1')) { + } else if (lsame_(norm, "O") || *(unsigned char *) + norm == '1') { /* Find norm1(A). */ @@ -14065,8 +14049,7 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, value = max(d__1,d__2); /* L80: */ } - } else if ((lsame_(norm, "F")) || (lsame_(norm, - "E"))) { + } else if (lsame_(norm, "F") || lsame_(norm, "E")) { /* Find normF(A). */ @@ -14176,7 +14159,7 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, --rwork; /* Function Body */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } @@ -14331,7 +14314,7 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, /* Form H * C */ - if ((tau->r != 0.) || (tau->i != 0.)) { + if (tau->r != 0. || tau->i != 0.) { /* w := C' * v */ @@ -14348,7 +14331,7 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, /* Form C * H */ - if ((tau->r != 0.) || (tau->i != 0.)) { + if (tau->r != 0. || tau->i != 0.) { /* w := C * v */ @@ -14495,7 +14478,7 @@ doublereal zlanhs_(char *norm, integer *n, doublecomplex *a, integer *lda, work -= work_offset; /* Function Body */ - if ((*m <= 0) || (*n <= 0)) { + if (*m <= 0 || *n <= 0) { return 0; } @@ -16771,7 +16754,7 @@ L230: i__2 = j + c_dim1; z__2.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__2.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__3.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__3.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__1.r = z__2.r + z__3.r, z__1.i = z__2.i + z__3.i; @@ -16782,8 +16765,8 @@ L230: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -16815,7 +16798,7 @@ L250: i__2 = j + c_dim1; z__3.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__3.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__4.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__4.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__2.r = z__3.r + z__4.r, z__2.i = z__3.i + z__4.i; @@ -16830,8 +16813,8 @@ L250: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -16874,7 +16857,7 @@ L270: i__2 = j + c_dim1; z__4.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__4.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__5.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__5.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__3.r = z__4.r + z__5.r, z__3.i = z__4.i + z__5.i; @@ -16882,7 +16865,7 @@ L270: z__6.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__6.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__2.r = z__3.r + z__6.r, z__2.i = z__3.i + z__6.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__7.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__7.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__1.r = z__2.r + z__7.r, z__1.i = z__2.i + z__7.i; @@ -16893,8 +16876,8 @@ L270: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -16905,8 +16888,8 @@ L270: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -16948,7 +16931,7 @@ L290: i__2 = j + c_dim1; z__5.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__5.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__6.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__6.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__4.r = z__5.r + z__6.r, z__4.i = z__5.i + z__6.i; @@ -16956,7 +16939,7 @@ L290: z__7.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__7.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__3.r = z__4.r + z__7.r, z__3.i = z__4.i + z__7.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__8.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__8.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__2.r = z__3.r + z__8.r, z__2.i = z__3.i + z__8.i; @@ -16971,8 +16954,8 @@ L290: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -16983,8 +16966,8 @@ L290: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17037,7 +17020,7 @@ L310: i__2 = j + c_dim1; z__6.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__6.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__7.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__7.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__5.r = z__6.r + z__7.r, z__5.i = z__6.i + z__7.i; @@ -17045,7 +17028,7 @@ L310: z__8.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__8.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__4.r = z__5.r + z__8.r, z__4.i = z__5.i + z__8.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__9.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__9.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__3.r = z__4.r + z__9.r, z__3.i = z__4.i + z__9.i; @@ -17064,8 +17047,8 @@ L310: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17076,8 +17059,8 @@ L310: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17141,7 +17124,7 @@ L330: i__2 = j + c_dim1; z__7.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__7.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__8.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__8.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__6.r = z__7.r + z__8.r, z__6.i = z__7.i + z__8.i; @@ -17149,7 +17132,7 @@ L330: z__9.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__9.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__5.r = z__6.r + z__9.r, z__5.i = z__6.i + z__9.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__10.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__10.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__4.r = z__5.r + z__10.r, z__4.i = z__5.i + z__10.i; @@ -17172,8 +17155,8 @@ L330: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17184,8 +17167,8 @@ L330: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17260,7 +17243,7 @@ L350: i__2 = j + c_dim1; z__8.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__8.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__9.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__9.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__7.r = z__8.r + z__9.r, z__7.i = z__8.i + z__9.i; @@ -17268,7 +17251,7 @@ L350: z__10.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__10.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__6.r = z__7.r + z__10.r, z__6.i = z__7.i + z__10.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__11.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__11.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__5.r = z__6.r + z__11.r, z__5.i = z__6.i + z__11.i; @@ -17284,7 +17267,7 @@ L350: z__14.r = v7.r * c__[i__8].r - v7.i * c__[i__8].i, z__14.i = v7.r * c__[i__8].i + v7.i * c__[i__8].r; z__2.r = z__3.r + z__14.r, z__2.i = z__3.i + z__14.i; - i__9 = j + ((c_dim1) << (3)); + i__9 = j + (c_dim1 << 3); z__15.r = v8.r * c__[i__9].r - v8.i * c__[i__9].i, z__15.i = v8.r * c__[i__9].i + v8.i * c__[i__9].r; z__1.r = z__2.r + z__15.r, z__1.i = z__2.i + z__15.i; @@ -17295,8 +17278,8 @@ L350: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17307,8 +17290,8 @@ L350: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17331,8 +17314,8 @@ L350: sum.i * t7.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (3)); - i__3 = j + ((c_dim1) << (3)); + i__2 = j + (c_dim1 << 3); + i__3 = j + (c_dim1 << 3); z__2.r = sum.r * t8.r - sum.i * t8.i, z__2.i = sum.r * t8.i + sum.i * t8.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17394,7 +17377,7 @@ L370: i__2 = j + c_dim1; z__9.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__9.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__10.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__10.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__8.r = z__9.r + z__10.r, z__8.i = z__9.i + z__10.i; @@ -17402,7 +17385,7 @@ L370: z__11.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__11.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__7.r = z__8.r + z__11.r, z__7.i = z__8.i + z__11.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__12.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__12.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__6.r = z__7.r + z__12.r, z__6.i = z__7.i + z__12.i; @@ -17418,7 +17401,7 @@ L370: z__15.r = v7.r * c__[i__8].r - v7.i * c__[i__8].i, z__15.i = v7.r * c__[i__8].i + v7.i * c__[i__8].r; z__3.r = z__4.r + z__15.r, z__3.i = z__4.i + z__15.i; - i__9 = j + ((c_dim1) << (3)); + i__9 = j + (c_dim1 << 3); z__16.r = v8.r * c__[i__9].r - v8.i * c__[i__9].i, z__16.i = v8.r * c__[i__9].i + v8.i * c__[i__9].r; z__2.r = z__3.r + z__16.r, z__2.i = z__3.i + z__16.i; @@ -17433,8 +17416,8 @@ L370: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17445,8 +17428,8 @@ L370: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17469,8 +17452,8 @@ L370: sum.i * t7.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (3)); - i__3 = j + ((c_dim1) << (3)); + i__2 = j + (c_dim1 << 3); + i__3 = j + (c_dim1 << 3); z__2.r = sum.r * t8.r - sum.i * t8.i, z__2.i = sum.r * t8.i + sum.i * t8.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17543,7 +17526,7 @@ L390: i__2 = j + c_dim1; z__10.r = v1.r * c__[i__2].r - v1.i * c__[i__2].i, z__10.i = v1.r * c__[i__2].i + v1.i * c__[i__2].r; - i__3 = j + ((c_dim1) << (1)); + i__3 = j + (c_dim1 << 1); z__11.r = v2.r * c__[i__3].r - v2.i * c__[i__3].i, z__11.i = v2.r * c__[i__3].i + v2.i * c__[i__3].r; z__9.r = z__10.r + z__11.r, z__9.i = z__10.i + z__11.i; @@ -17551,7 +17534,7 @@ L390: z__12.r = v3.r * c__[i__4].r - v3.i * c__[i__4].i, z__12.i = v3.r * c__[i__4].i + v3.i * c__[i__4].r; z__8.r = z__9.r + z__12.r, z__8.i = z__9.i + z__12.i; - i__5 = j + ((c_dim1) << (2)); + i__5 = j + (c_dim1 << 2); z__13.r = v4.r * c__[i__5].r - v4.i * c__[i__5].i, z__13.i = v4.r * c__[i__5].i + v4.i * c__[i__5].r; z__7.r = z__8.r + z__13.r, z__7.i = z__8.i + z__13.i; @@ -17567,7 +17550,7 @@ L390: z__16.r = v7.r * c__[i__8].r - v7.i * c__[i__8].i, z__16.i = v7.r * c__[i__8].i + v7.i * c__[i__8].r; z__4.r = z__5.r + z__16.r, z__4.i = z__5.i + z__16.i; - i__9 = j + ((c_dim1) << (3)); + i__9 = j + (c_dim1 << 3); z__17.r = v8.r * c__[i__9].r - v8.i * c__[i__9].i, z__17.i = v8.r * c__[i__9].i + v8.i * c__[i__9].r; z__3.r = z__4.r + z__17.r, z__3.i = z__4.i + z__17.i; @@ -17586,8 +17569,8 @@ L390: sum.i * t1.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (1)); - i__3 = j + ((c_dim1) << (1)); + i__2 = j + (c_dim1 << 1); + i__3 = j + (c_dim1 << 1); z__2.r = sum.r * t2.r - sum.i * t2.i, z__2.i = sum.r * t2.i + sum.i * t2.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17598,8 +17581,8 @@ L390: sum.i * t3.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (2)); - i__3 = j + ((c_dim1) << (2)); + i__2 = j + (c_dim1 << 2); + i__3 = j + (c_dim1 << 2); z__2.r = sum.r * t4.r - sum.i * t4.i, z__2.i = sum.r * t4.i + sum.i * t4.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17622,8 +17605,8 @@ L390: sum.i * t7.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; c__[i__2].r = z__1.r, c__[i__2].i = z__1.i; - i__2 = j + ((c_dim1) << (3)); - i__3 = j + ((c_dim1) << (3)); + i__2 = j + (c_dim1 << 3); + i__3 = j + (c_dim1 << 3); z__2.r = sum.r * t8.r - sum.i * t8.i, z__2.i = sum.r * t8.i + sum.i * t8.r; z__1.r = c__[i__3].r - z__2.r, z__1.i = c__[i__3].i - z__2.i; @@ -17777,25 +17760,23 @@ L410: *info = -4; } else if (*m < 0) { *info = -6; - } else if (((*n < 0) || (itype == 4 && *n != *m)) || (itype == 5 && *n != - *m)) { + } else if (*n < 0 || itype == 4 && *n != *m || itype == 5 && *n != *m) { *info = -7; } else if (itype <= 3 && *lda < max(1,*m)) { *info = -9; } else if (itype >= 4) { /* Computing MAX */ i__1 = *m - 1; - if ((*kl < 0) || (*kl > max(i__1,0))) { + if (*kl < 0 || *kl > max(i__1,0)) { *info = -2; } else /* if(complicated condition) */ { /* Computing MAX */ i__1 = *n - 1; - if (((*ku < 0) || (*ku > max(i__1,0))) || (((itype == 4) || ( - itype == 5)) && *kl != *ku)) { + if (*ku < 0 || *ku > max(i__1,0) || (itype == 4 || itype == 5) && + *kl != *ku) { *info = -3; - } else if (((itype == 4 && *lda < *kl + 1) || (itype == 5 && *lda - < *ku + 1)) || (itype == 6 && *lda < ((*kl) << (1)) + *ku - + 1)) { + } else if (itype == 4 && *lda < *kl + 1 || itype == 5 && *lda < * + ku + 1 || itype == 6 && *lda < (*kl << 1) + *ku + 1) { *info = -9; } } @@ -17809,7 +17790,7 @@ L410: /* Quick return if possible */ - if ((*n == 0) || (*m == 0)) { + if (*n == 0 || *m == 0) { return 0; } @@ -17955,7 +17936,7 @@ L10: k1 = *kl + *ku + 2; k2 = *kl + 1; - k3 = ((*kl) << (1)) + *ku + 1; + k3 = (*kl << 1) + *ku + 1; k4 = *kl + *ku + 1 + *m; i__1 = *n; for (j = 1; j <= i__1; ++j) { @@ -18257,13 +18238,13 @@ L10: /* Function Body */ info = 0; - if (! ((lsame_(side, "L")) || (lsame_(side, "R")))) { + if (! (lsame_(side, "L") || lsame_(side, "R"))) { info = 1; - } else if (! (((lsame_(pivot, "V")) || (lsame_( - pivot, "T"))) || (lsame_(pivot, "B")))) { + } else if (! (lsame_(pivot, "V") || lsame_(pivot, + "T") || lsame_(pivot, "B"))) { info = 2; - } else if (! ((lsame_(direct, "F")) || (lsame_( - direct, "B")))) { + } else if (! (lsame_(direct, "F") || lsame_(direct, + "B"))) { info = 3; } else if (*m < 0) { info = 4; @@ -18279,7 +18260,7 @@ L10: /* Quick return if possible */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } if (lsame_(side, "L")) { @@ -18292,7 +18273,7 @@ L10: for (j = 1; j <= i__1; ++j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = j + 1 + i__ * a_dim1; @@ -18322,7 +18303,7 @@ L10: for (j = *m - 1; j >= 1; --j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = j + 1 + i__ * a_dim1; @@ -18355,7 +18336,7 @@ L10: for (j = 2; j <= i__1; ++j) { ctemp = c__[j - 1]; stemp = s[j - 1]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = j + i__ * a_dim1; @@ -18385,7 +18366,7 @@ L10: for (j = *m; j >= 2; --j) { ctemp = c__[j - 1]; stemp = s[j - 1]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = j + i__ * a_dim1; @@ -18418,7 +18399,7 @@ L10: for (j = 1; j <= i__1; ++j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *n; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = j + i__ * a_dim1; @@ -18448,7 +18429,7 @@ L10: for (j = *m - 1; j >= 1; --j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = j + i__ * a_dim1; @@ -18486,7 +18467,7 @@ L10: for (j = 1; j <= i__1; ++j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = i__ + (j + 1) * a_dim1; @@ -18516,7 +18497,7 @@ L10: for (j = *n - 1; j >= 1; --j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + (j + 1) * a_dim1; @@ -18549,7 +18530,7 @@ L10: for (j = 2; j <= i__1; ++j) { ctemp = c__[j - 1]; stemp = s[j - 1]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = i__ + j * a_dim1; @@ -18579,7 +18560,7 @@ L10: for (j = *n; j >= 2; --j) { ctemp = c__[j - 1]; stemp = s[j - 1]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + j * a_dim1; @@ -18612,7 +18593,7 @@ L10: for (j = 1; j <= i__1; ++j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__2 = *m; for (i__ = 1; i__ <= i__2; ++i__) { i__3 = i__ + j * a_dim1; @@ -18642,7 +18623,7 @@ L10: for (j = *n - 1; j >= 1; --j) { ctemp = c__[j]; stemp = s[j]; - if ((ctemp != 1.) || (stemp != 0.)) { + if (ctemp != 1. || stemp != 0.) { i__1 = *m; for (i__ = 1; i__ <= i__1; ++i__) { i__2 = i__ + j * a_dim1; @@ -18882,7 +18863,7 @@ L10: return 0; } - n32 = (*n / 32) << (5); + n32 = *n / 32 << 5; if (n32 != 0) { i__1 = n32; for (j = 1; j <= i__1; j += 32) { @@ -20504,17 +20485,20 @@ L210: } /* zlatrs_ */ -/* Subroutine */ int zlauu2_(char *uplo, integer *n, doublecomplex *a, +/* Subroutine */ int zpotf2_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *info) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3; doublereal d__1; - doublecomplex z__1; + doublecomplex z__1, z__2; + + /* Builtin functions */ + double sqrt(doublereal); /* Local variables */ - static integer i__; - static doublereal aii; + static integer j; + static doublereal ajj; extern logical lsame_(char *, char *); extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, doublecomplex *, integer *, doublecomplex *, integer *); @@ -20528,7 +20512,7 @@ L210: /* - -- LAPACK auxiliary routine (version 3.0) -- + -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994 @@ -20537,35 +20521,39 @@ L210: Purpose ======= - ZLAUU2 computes the product U * U' or L' * L, where the triangular - factor U or L is stored in the upper or lower triangular part of - the array A. + ZPOTF2 computes the Cholesky factorization of a complex Hermitian + positive definite matrix A. - If UPLO = 'U' or 'u' then the upper triangle of the result is stored, - overwriting the factor U in A. - If UPLO = 'L' or 'l' then the lower triangle of the result is stored, - overwriting the factor L in A. + The factorization has the form + A = U' * U , if UPLO = 'U', or + A = L * L', if UPLO = 'L', + where U is an upper triangular matrix and L is lower triangular. - This is the unblocked form of the algorithm, calling Level 2 BLAS. + This is the unblocked version of the algorithm, calling Level 2 BLAS. Arguments ========= UPLO (input) CHARACTER*1 - Specifies whether the triangular factor stored in the array A - is upper or lower triangular: + Specifies whether the upper or lower triangular part of the + Hermitian matrix A is stored. = 'U': Upper triangular = 'L': Lower triangular N (input) INTEGER - The order of the triangular factor U or L. N >= 0. + The order of the matrix A. N >= 0. A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the triangular factor U or L. - On exit, if UPLO = 'U', the upper triangle of A is - overwritten with the upper triangle of the product U * U'; - if UPLO = 'L', the lower triangle of A is overwritten with - the lower triangle of the product L' * L. + On entry, the Hermitian matrix A. If UPLO = 'U', the leading + n by n upper triangular part of A contains the upper + triangular part of the matrix A, and the strictly lower + triangular part of A is not referenced. If UPLO = 'L', the + leading n by n lower triangular part of A contains the lower + triangular part of the matrix A, and the strictly upper + triangular part of A is not referenced. + + On exit, if INFO = 0, the factor U or L from the Cholesky + factorization A = U'*U or A = L*L'. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). @@ -20573,6 +20561,9 @@ L210: INFO (output) INTEGER = 0: successful exit < 0: if INFO = -k, the k-th argument had an illegal value + > 0: if INFO = k, the leading minor of order k is not + positive definite, and the factorization could not be + completed. ===================================================================== @@ -20597,7 +20588,7 @@ L210: } if (*info != 0) { i__1 = -(*info); - xerbla_("ZLAUU2", &i__1); + xerbla_("ZPOTF2", &i__1); return 0; } @@ -20609,81 +20600,113 @@ L210: if (upper) { -/* Compute the product U * U'. */ +/* Compute the Cholesky factorization A = U'*U. */ i__1 = *n; - for (i__ = 1; i__ <= i__1; ++i__) { - i__2 = i__ + i__ * a_dim1; - aii = a[i__2].r; - if (i__ < *n) { - i__2 = i__ + i__ * a_dim1; - i__3 = *n - i__; - zdotc_(&z__1, &i__3, &a[i__ + (i__ + 1) * a_dim1], lda, &a[ - i__ + (i__ + 1) * a_dim1], lda); - d__1 = aii * aii + z__1.r; - a[i__2].r = d__1, a[i__2].i = 0.; - i__2 = *n - i__; - zlacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda); - i__2 = i__ - 1; - i__3 = *n - i__; - z__1.r = aii, z__1.i = 0.; - zgemv_("No transpose", &i__2, &i__3, &c_b60, &a[(i__ + 1) * - a_dim1 + 1], lda, &a[i__ + (i__ + 1) * a_dim1], lda, & - z__1, &a[i__ * a_dim1 + 1], &c__1); - i__2 = *n - i__; - zlacgv_(&i__2, &a[i__ + (i__ + 1) * a_dim1], lda); - } else { - zdscal_(&i__, &aii, &a[i__ * a_dim1 + 1], &c__1); + for (j = 1; j <= i__1; ++j) { + +/* Compute U(J,J) and test for non-positive-definiteness. */ + + i__2 = j + j * a_dim1; + d__1 = a[i__2].r; + i__3 = j - 1; + zdotc_(&z__2, &i__3, &a[j * a_dim1 + 1], &c__1, &a[j * a_dim1 + 1] + , &c__1); + z__1.r = d__1 - z__2.r, z__1.i = -z__2.i; + ajj = z__1.r; + if (ajj <= 0.) { + i__2 = j + j * a_dim1; + a[i__2].r = ajj, a[i__2].i = 0.; + goto L30; + } + ajj = sqrt(ajj); + i__2 = j + j * a_dim1; + a[i__2].r = ajj, a[i__2].i = 0.; + +/* Compute elements J+1:N of row J. */ + + if (j < *n) { + i__2 = j - 1; + zlacgv_(&i__2, &a[j * a_dim1 + 1], &c__1); + i__2 = j - 1; + i__3 = *n - j; + z__1.r = -1., z__1.i = -0.; + zgemv_("Transpose", &i__2, &i__3, &z__1, &a[(j + 1) * a_dim1 + + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b60, &a[j + ( + j + 1) * a_dim1], lda); + i__2 = j - 1; + zlacgv_(&i__2, &a[j * a_dim1 + 1], &c__1); + i__2 = *n - j; + d__1 = 1. / ajj; + zdscal_(&i__2, &d__1, &a[j + (j + 1) * a_dim1], lda); } /* L10: */ } - } else { -/* Compute the product L' * L. */ +/* Compute the Cholesky factorization A = L*L'. */ i__1 = *n; - for (i__ = 1; i__ <= i__1; ++i__) { - i__2 = i__ + i__ * a_dim1; - aii = a[i__2].r; - if (i__ < *n) { - i__2 = i__ + i__ * a_dim1; - i__3 = *n - i__; - zdotc_(&z__1, &i__3, &a[i__ + 1 + i__ * a_dim1], &c__1, &a[ - i__ + 1 + i__ * a_dim1], &c__1); - d__1 = aii * aii + z__1.r; - a[i__2].r = d__1, a[i__2].i = 0.; - i__2 = i__ - 1; - zlacgv_(&i__2, &a[i__ + a_dim1], lda); - i__2 = *n - i__; - i__3 = i__ - 1; - z__1.r = aii, z__1.i = 0.; - zgemv_("Conjugate transpose", &i__2, &i__3, &c_b60, &a[i__ + - 1 + a_dim1], lda, &a[i__ + 1 + i__ * a_dim1], &c__1, & - z__1, &a[i__ + a_dim1], lda); - i__2 = i__ - 1; - zlacgv_(&i__2, &a[i__ + a_dim1], lda); - } else { - zdscal_(&i__, &aii, &a[i__ + a_dim1], lda); + for (j = 1; j <= i__1; ++j) { + +/* Compute L(J,J) and test for non-positive-definiteness. */ + + i__2 = j + j * a_dim1; + d__1 = a[i__2].r; + i__3 = j - 1; + zdotc_(&z__2, &i__3, &a[j + a_dim1], lda, &a[j + a_dim1], lda); + z__1.r = d__1 - z__2.r, z__1.i = -z__2.i; + ajj = z__1.r; + if (ajj <= 0.) { + i__2 = j + j * a_dim1; + a[i__2].r = ajj, a[i__2].i = 0.; + goto L30; + } + ajj = sqrt(ajj); + i__2 = j + j * a_dim1; + a[i__2].r = ajj, a[i__2].i = 0.; + +/* Compute elements J+1:N of column J. */ + + if (j < *n) { + i__2 = j - 1; + zlacgv_(&i__2, &a[j + a_dim1], lda); + i__2 = *n - j; + i__3 = j - 1; + z__1.r = -1., z__1.i = -0.; + zgemv_("No transpose", &i__2, &i__3, &z__1, &a[j + 1 + a_dim1] + , lda, &a[j + a_dim1], lda, &c_b60, &a[j + 1 + j * + a_dim1], &c__1); + i__2 = j - 1; + zlacgv_(&i__2, &a[j + a_dim1], lda); + i__2 = *n - j; + d__1 = 1. / ajj; + zdscal_(&i__2, &d__1, &a[j + 1 + j * a_dim1], &c__1); } /* L20: */ } } + goto L40; +L30: + *info = j; + +L40: return 0; -/* End of ZLAUU2 */ +/* End of ZPOTF2 */ -} /* zlauu2_ */ +} /* zpotf2_ */ -/* Subroutine */ int zlauum_(char *uplo, integer *n, doublecomplex *a, +/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, integer *lda, integer *info) { /* System generated locals */ integer a_dim1, a_offset, i__1, i__2, i__3, i__4; + doublecomplex z__1; /* Local variables */ - static integer i__, ib, nb; + static integer j, jb, nb; extern logical lsame_(char *, char *); extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, @@ -20692,16 +20715,16 @@ L210: integer *, doublereal *, doublecomplex *, integer *, doublereal *, doublecomplex *, integer *); static logical upper; - extern /* Subroutine */ int ztrmm_(char *, char *, char *, char *, + extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *), - zlauu2_(char *, integer *, doublecomplex *, integer *, integer *), xerbla_(char *, integer *); + zpotf2_(char *, integer *, doublecomplex *, integer *, integer *), xerbla_(char *, integer *); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); /* - -- LAPACK auxiliary routine (version 3.0) -- + -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University September 30, 1994 @@ -20710,42 +20733,47 @@ L210: Purpose ======= - ZLAUUM computes the product U * U' or L' * L, where the triangular - factor U or L is stored in the upper or lower triangular part of - the array A. + ZPOTRF computes the Cholesky factorization of a complex Hermitian + positive definite matrix A. - If UPLO = 'U' or 'u' then the upper triangle of the result is stored, - overwriting the factor U in A. - If UPLO = 'L' or 'l' then the lower triangle of the result is stored, - overwriting the factor L in A. + The factorization has the form + A = U**H * U, if UPLO = 'U', or + A = L * L**H, if UPLO = 'L', + where U is an upper triangular matrix and L is lower triangular. - This is the blocked form of the algorithm, calling Level 3 BLAS. + This is the block version of the algorithm, calling Level 3 BLAS. Arguments ========= UPLO (input) CHARACTER*1 - Specifies whether the triangular factor stored in the array A - is upper or lower triangular: - = 'U': Upper triangular - = 'L': Lower triangular + = 'U': Upper triangle of A is stored; + = 'L': Lower triangle of A is stored. N (input) INTEGER - The order of the triangular factor U or L. N >= 0. + The order of the matrix A. N >= 0. A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the triangular factor U or L. - On exit, if UPLO = 'U', the upper triangle of A is - overwritten with the upper triangle of the product U * U'; - if UPLO = 'L', the lower triangle of A is overwritten with - the lower triangle of the product L' * L. + On entry, the Hermitian matrix A. If UPLO = 'U', the leading + N-by-N upper triangular part of A contains the upper + triangular part of the matrix A, and the strictly lower + triangular part of A is not referenced. If UPLO = 'L', the + leading N-by-N lower triangular part of A contains the lower + triangular part of the matrix A, and the strictly upper + triangular part of A is not referenced. + + On exit, if INFO = 0, the factor U or L from the Cholesky + factorization A = U**H*U or A = L*L**H. LDA (input) INTEGER The leading dimension of the array A. LDA >= max(1,N). INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -k, the k-th argument had an illegal value + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value + > 0: if INFO = i, the leading minor of order i is not + positive definite, and the factorization could not be + completed. ===================================================================== @@ -20770,7 +20798,7 @@ L210: } if (*info != 0) { i__1 = -(*info); - xerbla_("ZLAUUM", &i__1); + xerbla_("ZPOTRF", &i__1); return 0; } @@ -20782,164 +20810,285 @@ L210: /* Determine the block size for this environment. */ - nb = ilaenv_(&c__1, "ZLAUUM", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( + nb = ilaenv_(&c__1, "ZPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( ftnlen)1); + if (nb <= 1 || nb >= *n) { - if ((nb <= 1) || (nb >= *n)) { - -/* Use unblocked code */ +/* Use unblocked code. */ - zlauu2_(uplo, n, &a[a_offset], lda, info); + zpotf2_(uplo, n, &a[a_offset], lda, info); } else { -/* Use blocked code */ +/* Use blocked code. */ if (upper) { -/* Compute the product U * U'. */ +/* Compute the Cholesky factorization A = U'*U. */ i__1 = *n; i__2 = nb; - for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { -/* Computing MIN */ - i__3 = nb, i__4 = *n - i__ + 1; - ib = min(i__3,i__4); - i__3 = i__ - 1; - ztrmm_("Right", "Upper", "Conjugate transpose", "Non-unit", & - i__3, &ib, &c_b60, &a[i__ + i__ * a_dim1], lda, &a[ - i__ * a_dim1 + 1], lda); - zlauu2_("Upper", &ib, &a[i__ + i__ * a_dim1], lda, info); - if (i__ + ib <= *n) { - i__3 = i__ - 1; - i__4 = *n - i__ - ib + 1; - zgemm_("No transpose", "Conjugate transpose", &i__3, &ib, - &i__4, &c_b60, &a[(i__ + ib) * a_dim1 + 1], lda, & - a[i__ + (i__ + ib) * a_dim1], lda, &c_b60, &a[i__ - * a_dim1 + 1], lda); - i__3 = *n - i__ - ib + 1; - zherk_("Upper", "No transpose", &ib, &i__3, &c_b1015, &a[ - i__ + (i__ + ib) * a_dim1], lda, &c_b1015, &a[i__ - + i__ * a_dim1], lda); + for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { + +/* + Update and factorize the current diagonal block and test + for non-positive-definiteness. + + Computing MIN +*/ + i__3 = nb, i__4 = *n - j + 1; + jb = min(i__3,i__4); + i__3 = j - 1; + zherk_("Upper", "Conjugate transpose", &jb, &i__3, &c_b1294, & + a[j * a_dim1 + 1], lda, &c_b1015, &a[j + j * a_dim1], + lda); + zpotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info); + if (*info != 0) { + goto L30; } -/* L10: */ - } - } else { + if (j + jb <= *n) { -/* Compute the product L' * L. */ +/* Compute the current block row. */ - i__2 = *n; - i__1 = nb; - for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) { -/* Computing MIN */ - i__3 = nb, i__4 = *n - i__ + 1; - ib = min(i__3,i__4); - i__3 = i__ - 1; - ztrmm_("Left", "Lower", "Conjugate transpose", "Non-unit", & - ib, &i__3, &c_b60, &a[i__ + i__ * a_dim1], lda, &a[ - i__ + a_dim1], lda); - zlauu2_("Lower", &ib, &a[i__ + i__ * a_dim1], lda, info); - if (i__ + ib <= *n) { - i__3 = i__ - 1; - i__4 = *n - i__ - ib + 1; - zgemm_("Conjugate transpose", "No transpose", &ib, &i__3, - &i__4, &c_b60, &a[i__ + ib + i__ * a_dim1], lda, & - a[i__ + ib + a_dim1], lda, &c_b60, &a[i__ + + i__3 = *n - j - jb + 1; + i__4 = j - 1; + z__1.r = -1., z__1.i = -0.; + zgemm_("Conjugate transpose", "No transpose", &jb, &i__3, + &i__4, &z__1, &a[j * a_dim1 + 1], lda, &a[(j + jb) + * a_dim1 + 1], lda, &c_b60, &a[j + (j + jb) * a_dim1], lda); - i__3 = *n - i__ - ib + 1; - zherk_("Lower", "Conjugate transpose", &ib, &i__3, & - c_b1015, &a[i__ + ib + i__ * a_dim1], lda, & - c_b1015, &a[i__ + i__ * a_dim1], lda); + i__3 = *n - j - jb + 1; + ztrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", + &jb, &i__3, &c_b60, &a[j + j * a_dim1], lda, &a[ + j + (j + jb) * a_dim1], lda); } -/* L20: */ +/* L10: */ } - } - } - return 0; + } else { -/* End of ZLAUUM */ +/* Compute the Cholesky factorization A = L*L'. */ -} /* zlauum_ */ + i__2 = *n; + i__1 = nb; + for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { -/* Subroutine */ int zpotf2_(char *uplo, integer *n, doublecomplex *a, - integer *lda, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2, i__3; - doublereal d__1; - doublecomplex z__1, z__2; +/* + Update and factorize the current diagonal block and test + for non-positive-definiteness. - /* Builtin functions */ + Computing MIN +*/ + i__3 = nb, i__4 = *n - j + 1; + jb = min(i__3,i__4); + i__3 = j - 1; + zherk_("Lower", "No transpose", &jb, &i__3, &c_b1294, &a[j + + a_dim1], lda, &c_b1015, &a[j + j * a_dim1], lda); + zpotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info); + if (*info != 0) { + goto L30; + } + if (j + jb <= *n) { + +/* Compute the current block column. */ + + i__3 = *n - j - jb + 1; + i__4 = j - 1; + z__1.r = -1., z__1.i = -0.; + zgemm_("No transpose", "Conjugate transpose", &i__3, &jb, + &i__4, &z__1, &a[j + jb + a_dim1], lda, &a[j + + a_dim1], lda, &c_b60, &a[j + jb + j * a_dim1], + lda); + i__3 = *n - j - jb + 1; + ztrsm_("Right", "Lower", "Conjugate transpose", "Non-unit" + , &i__3, &jb, &c_b60, &a[j + j * a_dim1], lda, &a[ + j + jb + j * a_dim1], lda); + } +/* L20: */ + } + } + } + goto L40; + +L30: + *info = *info + j - 1; + +L40: + return 0; + +/* End of ZPOTRF */ + +} /* zpotrf_ */ + +/* Subroutine */ int zstedc_(char *compz, integer *n, doublereal *d__, + doublereal *e, doublecomplex *z__, integer *ldz, doublecomplex *work, + integer *lwork, doublereal *rwork, integer *lrwork, integer *iwork, + integer *liwork, integer *info) +{ + /* System generated locals */ + integer z_dim1, z_offset, i__1, i__2, i__3, i__4; + doublereal d__1, d__2; + + /* Builtin functions */ + double log(doublereal); + integer pow_ii(integer *, integer *); double sqrt(doublereal); /* Local variables */ - static integer j; - static doublereal ajj; + static integer i__, j, k, m; + static doublereal p; + static integer ii, ll, end, lgn; + static doublereal eps, tiny; extern logical lsame_(char *, char *); - extern /* Double Complex */ VOID zdotc_(doublecomplex *, integer *, - doublecomplex *, integer *, doublecomplex *, integer *); - extern /* Subroutine */ int zgemv_(char *, integer *, integer *, - doublecomplex *, doublecomplex *, integer *, doublecomplex *, - integer *, doublecomplex *, doublecomplex *, integer *); - static logical upper; - extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_( - integer *, doublereal *, doublecomplex *, integer *), zlacgv_( + static integer lwmin, start; + extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, + doublecomplex *, integer *), zlaed0_(integer *, integer *, + doublereal *, doublereal *, doublecomplex *, integer *, + doublecomplex *, integer *, doublereal *, integer *, integer *); + + extern /* Subroutine */ int dlascl_(char *, integer *, integer *, + doublereal *, doublereal *, integer *, integer *, doublereal *, + integer *, integer *), dstedc_(char *, integer *, + doublereal *, doublereal *, doublereal *, integer *, doublereal *, + integer *, integer *, integer *, integer *), dlaset_( + char *, integer *, integer *, doublereal *, doublereal *, + doublereal *, integer *), xerbla_(char *, integer *); + extern integer ilaenv_(integer *, char *, char *, integer *, integer *, + integer *, integer *, ftnlen, ftnlen); + extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *); + extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *, + integer *), zlacrm_(integer *, integer *, doublecomplex *, + integer *, doublereal *, integer *, doublecomplex *, integer *, + doublereal *); + static integer liwmin, icompz; + extern /* Subroutine */ int dsteqr_(char *, integer *, doublereal *, + doublereal *, doublereal *, integer *, doublereal *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *); + static doublereal orgnrm; + static integer lrwmin; + static logical lquery; + static integer smlsiz; + extern /* Subroutine */ int zsteqr_(char *, integer *, doublereal *, + doublereal *, doublecomplex *, integer *, doublereal *, integer *); /* -- LAPACK routine (version 3.0) -- Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 + June 30, 1999 Purpose ======= - ZPOTF2 computes the Cholesky factorization of a complex Hermitian - positive definite matrix A. - - The factorization has the form - A = U' * U , if UPLO = 'U', or - A = L * L', if UPLO = 'L', - where U is an upper triangular matrix and L is lower triangular. + ZSTEDC computes all eigenvalues and, optionally, eigenvectors of a + symmetric tridiagonal matrix using the divide and conquer method. + The eigenvectors of a full or band complex Hermitian matrix can also + be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this + matrix to tridiagonal form. - This is the unblocked version of the algorithm, calling Level 2 BLAS. + This code makes very mild assumptions about floating point + arithmetic. It will work on machines with a guard digit in + add/subtract, or on those binary machines without guard digits + which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. + It could conceivably fail on hexadecimal or decimal machines + without guard digits, but we know of none. See DLAED3 for details. Arguments ========= - UPLO (input) CHARACTER*1 - Specifies whether the upper or lower triangular part of the - Hermitian matrix A is stored. - = 'U': Upper triangular - = 'L': Lower triangular + COMPZ (input) CHARACTER*1 + = 'N': Compute eigenvalues only. + = 'I': Compute eigenvectors of tridiagonal matrix also. + = 'V': Compute eigenvectors of original Hermitian matrix + also. On entry, Z contains the unitary matrix used + to reduce the original matrix to tridiagonal form. N (input) INTEGER - The order of the matrix A. N >= 0. + The dimension of the symmetric tridiagonal matrix. N >= 0. - A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the Hermitian matrix A. If UPLO = 'U', the leading - n by n upper triangular part of A contains the upper - triangular part of the matrix A, and the strictly lower - triangular part of A is not referenced. If UPLO = 'L', the - leading n by n lower triangular part of A contains the lower - triangular part of the matrix A, and the strictly upper - triangular part of A is not referenced. + D (input/output) DOUBLE PRECISION array, dimension (N) + On entry, the diagonal elements of the tridiagonal matrix. + On exit, if INFO = 0, the eigenvalues in ascending order. - On exit, if INFO = 0, the factor U or L from the Cholesky - factorization A = U'*U or A = L*L'. + E (input/output) DOUBLE PRECISION array, dimension (N-1) + On entry, the subdiagonal elements of the tridiagonal matrix. + On exit, E has been destroyed. - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). + Z (input/output) COMPLEX*16 array, dimension (LDZ,N) + On entry, if COMPZ = 'V', then Z contains the unitary + matrix used in the reduction to tridiagonal form. + On exit, if INFO = 0, then if COMPZ = 'V', Z contains the + orthonormal eigenvectors of the original Hermitian matrix, + and if COMPZ = 'I', Z contains the orthonormal eigenvectors + of the symmetric tridiagonal matrix. + If COMPZ = 'N', then Z is not referenced. + + LDZ (input) INTEGER + The leading dimension of the array Z. LDZ >= 1. + If eigenvectors are desired, then LDZ >= max(1,N). + + WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) + On exit, if INFO = 0, WORK(1) returns the optimal LWORK. + + LWORK (input) INTEGER + The dimension of the array WORK. + If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1. + If COMPZ = 'V' and N > 1, LWORK must be at least N*N. + + If LWORK = -1, then a workspace query is assumed; the routine + only calculates the optimal size of the WORK array, returns + this value as the first entry of the WORK array, and no error + message related to LWORK is issued by XERBLA. + + RWORK (workspace/output) DOUBLE PRECISION array, + dimension (LRWORK) + On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. + + LRWORK (input) INTEGER + The dimension of the array RWORK. + If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. + If COMPZ = 'V' and N > 1, LRWORK must be at least + 1 + 3*N + 2*N*lg N + 3*N**2 , + where lg( N ) = smallest integer k such + that 2**k >= N. + If COMPZ = 'I' and N > 1, LRWORK must be at least + 1 + 4*N + 2*N**2 . + + If LRWORK = -1, then a workspace query is assumed; the + routine only calculates the optimal size of the RWORK array, + returns this value as the first entry of the RWORK array, and + no error message related to LRWORK is issued by XERBLA. + + IWORK (workspace/output) INTEGER array, dimension (LIWORK) + On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. + + LIWORK (input) INTEGER + The dimension of the array IWORK. + If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. + If COMPZ = 'V' or N > 1, LIWORK must be at least + 6 + 6*N + 5*N*lg N. + If COMPZ = 'I' or N > 1, LIWORK must be at least + 3 + 5*N . + + If LIWORK = -1, then a workspace query is assumed; the + routine only calculates the optimal size of the IWORK array, + returns this value as the first entry of the IWORK array, and + no error message related to LIWORK is issued by XERBLA. INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -k, the k-th argument had an illegal value - > 0: if INFO = k, the leading minor of order k is not - positive definite, and the factorization could not be - completed. + = 0: successful exit. + < 0: if INFO = -i, the i-th argument had an illegal value. + > 0: The algorithm failed to compute an eigenvalue while + working on the submatrix lying in rows and columns + INFO/(N+1) through mod(INFO,N+1). + + Further Details + =============== + + Based on contributions by + Jeff Rutter, Computer Science Division, University of California + at Berkeley, USA ===================================================================== @@ -20948,208 +21097,403 @@ L210: */ /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; + --d__; + --e; + z_dim1 = *ldz; + z_offset = 1 + z_dim1; + z__ -= z_offset; + --work; + --rwork; + --iwork; /* Function Body */ *info = 0; - upper = lsame_(uplo, "U"); - if (! upper && ! lsame_(uplo, "L")) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*lda < max(1,*n)) { - *info = -4; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZPOTF2", &i__1); - return 0; - } - -/* Quick return if possible */ + lquery = *lwork == -1 || *lrwork == -1 || *liwork == -1; - if (*n == 0) { - return 0; - } - - if (upper) { + if (lsame_(compz, "N")) { + icompz = 0; + } else if (lsame_(compz, "V")) { + icompz = 1; + } else if (lsame_(compz, "I")) { + icompz = 2; + } else { + icompz = -1; + } + if (*n <= 1 || icompz <= 0) { + lwmin = 1; + liwmin = 1; + lrwmin = 1; + } else { + lgn = (integer) (log((doublereal) (*n)) / log(2.)); + if (pow_ii(&c__2, &lgn) < *n) { + ++lgn; + } + if (pow_ii(&c__2, &lgn) < *n) { + ++lgn; + } + if (icompz == 1) { + lwmin = *n * *n; +/* Computing 2nd power */ + i__1 = *n; + lrwmin = *n * 3 + 1 + (*n << 1) * lgn + i__1 * i__1 * 3; + liwmin = *n * 6 + 6 + *n * 5 * lgn; + } else if (icompz == 2) { + lwmin = 1; +/* Computing 2nd power */ + i__1 = *n; + lrwmin = (*n << 2) + 1 + (i__1 * i__1 << 1); + liwmin = *n * 5 + 3; + } + } + if (icompz < 0) { + *info = -1; + } else if (*n < 0) { + *info = -2; + } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) { + *info = -6; + } else if (*lwork < lwmin && ! lquery) { + *info = -8; + } else if (*lrwork < lrwmin && ! lquery) { + *info = -10; + } else if (*liwork < liwmin && ! lquery) { + *info = -12; + } -/* Compute the Cholesky factorization A = U'*U. */ + if (*info == 0) { + work[1].r = (doublereal) lwmin, work[1].i = 0.; + rwork[1] = (doublereal) lrwmin; + iwork[1] = liwmin; + } - i__1 = *n; - for (j = 1; j <= i__1; ++j) { + if (*info != 0) { + i__1 = -(*info); + xerbla_("ZSTEDC", &i__1); + return 0; + } else if (lquery) { + return 0; + } -/* Compute U(J,J) and test for non-positive-definiteness. */ +/* Quick return if possible */ - i__2 = j + j * a_dim1; - d__1 = a[i__2].r; - i__3 = j - 1; - zdotc_(&z__2, &i__3, &a[j * a_dim1 + 1], &c__1, &a[j * a_dim1 + 1] - , &c__1); - z__1.r = d__1 - z__2.r, z__1.i = -z__2.i; - ajj = z__1.r; - if (ajj <= 0.) { - i__2 = j + j * a_dim1; - a[i__2].r = ajj, a[i__2].i = 0.; - goto L30; - } - ajj = sqrt(ajj); - i__2 = j + j * a_dim1; - a[i__2].r = ajj, a[i__2].i = 0.; + if (*n == 0) { + return 0; + } + if (*n == 1) { + if (icompz != 0) { + i__1 = z_dim1 + 1; + z__[i__1].r = 1., z__[i__1].i = 0.; + } + return 0; + } -/* Compute elements J+1:N of row J. */ + smlsiz = ilaenv_(&c__9, "ZSTEDC", " ", &c__0, &c__0, &c__0, &c__0, ( + ftnlen)6, (ftnlen)1); - if (j < *n) { - i__2 = j - 1; - zlacgv_(&i__2, &a[j * a_dim1 + 1], &c__1); - i__2 = j - 1; - i__3 = *n - j; - z__1.r = -1., z__1.i = -0.; - zgemv_("Transpose", &i__2, &i__3, &z__1, &a[(j + 1) * a_dim1 - + 1], lda, &a[j * a_dim1 + 1], &c__1, &c_b60, &a[j + ( - j + 1) * a_dim1], lda); - i__2 = j - 1; - zlacgv_(&i__2, &a[j * a_dim1 + 1], &c__1); - i__2 = *n - j; - d__1 = 1. / ajj; - zdscal_(&i__2, &d__1, &a[j + (j + 1) * a_dim1], lda); - } -/* L10: */ - } - } else { +/* + If the following conditional clause is removed, then the routine + will use the Divide and Conquer routine to compute only the + eigenvalues, which requires (3N + 3N**2) real workspace and + (2 + 5N + 2N lg(N)) integer workspace. + Since on many architectures DSTERF is much faster than any other + algorithm for finding eigenvalues only, it is used here + as the default. -/* Compute the Cholesky factorization A = L*L'. */ + If COMPZ = 'N', use DSTERF to compute the eigenvalues. +*/ - i__1 = *n; - for (j = 1; j <= i__1; ++j) { + if (icompz == 0) { + dsterf_(n, &d__[1], &e[1], info); + return 0; + } -/* Compute L(J,J) and test for non-positive-definiteness. */ +/* + If N is smaller than the minimum divide size (SMLSIZ+1), then + solve the problem with another solver. +*/ - i__2 = j + j * a_dim1; - d__1 = a[i__2].r; - i__3 = j - 1; - zdotc_(&z__2, &i__3, &a[j + a_dim1], lda, &a[j + a_dim1], lda); - z__1.r = d__1 - z__2.r, z__1.i = -z__2.i; - ajj = z__1.r; - if (ajj <= 0.) { - i__2 = j + j * a_dim1; - a[i__2].r = ajj, a[i__2].i = 0.; - goto L30; - } - ajj = sqrt(ajj); - i__2 = j + j * a_dim1; - a[i__2].r = ajj, a[i__2].i = 0.; + if (*n <= smlsiz) { + if (icompz == 0) { + dsterf_(n, &d__[1], &e[1], info); + return 0; + } else if (icompz == 2) { + zsteqr_("I", n, &d__[1], &e[1], &z__[z_offset], ldz, &rwork[1], + info); + return 0; + } else { + zsteqr_("V", n, &d__[1], &e[1], &z__[z_offset], ldz, &rwork[1], + info); + return 0; + } + } -/* Compute elements J+1:N of column J. */ +/* If COMPZ = 'I', we simply call DSTEDC instead. */ - if (j < *n) { - i__2 = j - 1; - zlacgv_(&i__2, &a[j + a_dim1], lda); - i__2 = *n - j; - i__3 = j - 1; - z__1.r = -1., z__1.i = -0.; - zgemv_("No transpose", &i__2, &i__3, &z__1, &a[j + 1 + a_dim1] - , lda, &a[j + a_dim1], lda, &c_b60, &a[j + 1 + j * - a_dim1], &c__1); - i__2 = j - 1; - zlacgv_(&i__2, &a[j + a_dim1], lda); - i__2 = *n - j; - d__1 = 1. / ajj; - zdscal_(&i__2, &d__1, &a[j + 1 + j * a_dim1], &c__1); + if (icompz == 2) { + dlaset_("Full", n, n, &c_b324, &c_b1015, &rwork[1], n); + ll = *n * *n + 1; + i__1 = *lrwork - ll + 1; + dstedc_("I", n, &d__[1], &e[1], &rwork[1], n, &rwork[ll], &i__1, & + iwork[1], liwork, info); + i__1 = *n; + for (j = 1; j <= i__1; ++j) { + i__2 = *n; + for (i__ = 1; i__ <= i__2; ++i__) { + i__3 = i__ + j * z_dim1; + i__4 = (j - 1) * *n + i__; + z__[i__3].r = rwork[i__4], z__[i__3].i = 0.; +/* L10: */ } /* L20: */ } + return 0; } - goto L40; -L30: - *info = j; +/* + From now on, only option left to be handled is COMPZ = 'V', + i.e. ICOMPZ = 1. -L40: - return 0; + Scale. +*/ -/* End of ZPOTF2 */ + orgnrm = dlanst_("M", n, &d__[1], &e[1]); + if (orgnrm == 0.) { + return 0; + } -} /* zpotf2_ */ + eps = EPSILON; -/* Subroutine */ int zpotrf_(char *uplo, integer *n, doublecomplex *a, - integer *lda, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2, i__3, i__4; - doublecomplex z__1; + start = 1; - /* Local variables */ - static integer j, jb, nb; - extern logical lsame_(char *, char *); - extern /* Subroutine */ int zgemm_(char *, char *, integer *, integer *, - integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *, doublecomplex *, doublecomplex *, - integer *), zherk_(char *, char *, integer *, - integer *, doublereal *, doublecomplex *, integer *, doublereal *, - doublecomplex *, integer *); - static logical upper; - extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, - integer *, integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *), - zpotf2_(char *, integer *, doublecomplex *, integer *, integer *), xerbla_(char *, integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); +/* while ( START <= N ) */ +L30: + if (start <= *n) { /* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 + Let END be the position of the next subdiagonal entry such that + E( END ) <= TINY or END = N if no such subdiagonal exists. The + matrix identified by the elements between START and END + constitutes an independent sub-problem. +*/ + end = start; +L40: + if (end < *n) { + tiny = eps * sqrt((d__1 = d__[end], abs(d__1))) * sqrt((d__2 = + d__[end + 1], abs(d__2))); + if ((d__1 = e[end], abs(d__1)) > tiny) { + ++end; + goto L40; + } + } - Purpose - ======= +/* (Sub) Problem determined. Compute its size and solve it. */ - ZPOTRF computes the Cholesky factorization of a complex Hermitian - positive definite matrix A. + m = end - start + 1; + if (m > smlsiz) { + *info = smlsiz; - The factorization has the form - A = U**H * U, if UPLO = 'U', or - A = L * L**H, if UPLO = 'L', - where U is an upper triangular matrix and L is lower triangular. +/* Scale. */ - This is the block version of the algorithm, calling Level 3 BLAS. + orgnrm = dlanst_("M", &m, &d__[start], &e[start]); + dlascl_("G", &c__0, &c__0, &orgnrm, &c_b1015, &m, &c__1, &d__[ + start], &m, info); + i__1 = m - 1; + i__2 = m - 1; + dlascl_("G", &c__0, &c__0, &orgnrm, &c_b1015, &i__1, &c__1, &e[ + start], &i__2, info); + + zlaed0_(n, &m, &d__[start], &e[start], &z__[start * z_dim1 + 1], + ldz, &work[1], n, &rwork[1], &iwork[1], info); + if (*info > 0) { + *info = (*info / (m + 1) + start - 1) * (*n + 1) + *info % (m + + 1) + start - 1; + return 0; + } + +/* Scale back. */ + + dlascl_("G", &c__0, &c__0, &c_b1015, &orgnrm, &m, &c__1, &d__[ + start], &m, info); + + } else { + dsteqr_("I", &m, &d__[start], &e[start], &rwork[1], &m, &rwork[m * + m + 1], info); + zlacrm_(n, &m, &z__[start * z_dim1 + 1], ldz, &rwork[1], &m, & + work[1], n, &rwork[m * m + 1]); + zlacpy_("A", n, &m, &work[1], n, &z__[start * z_dim1 + 1], ldz); + if (*info > 0) { + *info = start * (*n + 1) + end; + return 0; + } + } + + start = end + 1; + goto L30; + } + +/* + endwhile + + If the problem split any number of times, then the eigenvalues + will not be properly ordered. Here we permute the eigenvalues + (and the associated eigenvectors) into ascending order. +*/ + + if (m != *n) { + +/* Use Selection Sort to minimize swaps of eigenvectors */ + + i__1 = *n; + for (ii = 2; ii <= i__1; ++ii) { + i__ = ii - 1; + k = i__; + p = d__[i__]; + i__2 = *n; + for (j = ii; j <= i__2; ++j) { + if (d__[j] < p) { + k = j; + p = d__[j]; + } +/* L50: */ + } + if (k != i__) { + d__[k] = d__[i__]; + d__[i__] = p; + zswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1], + &c__1); + } +/* L60: */ + } + } + + work[1].r = (doublereal) lwmin, work[1].i = 0.; + rwork[1] = (doublereal) lrwmin; + iwork[1] = liwmin; + + return 0; + +/* End of ZSTEDC */ + +} /* zstedc_ */ + +/* Subroutine */ int zsteqr_(char *compz, integer *n, doublereal *d__, + doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work, + integer *info) +{ + /* System generated locals */ + integer z_dim1, z_offset, i__1, i__2; + doublereal d__1, d__2; + + /* Builtin functions */ + double sqrt(doublereal), d_sign(doublereal *, doublereal *); + + /* Local variables */ + static doublereal b, c__, f, g; + static integer i__, j, k, l, m; + static doublereal p, r__, s; + static integer l1, ii, mm, lm1, mm1, nm1; + static doublereal rt1, rt2, eps; + static integer lsv; + static doublereal tst, eps2; + static integer lend, jtot; + extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal + *, doublereal *, doublereal *); + extern logical lsame_(char *, char *); + static doublereal anorm; + extern /* Subroutine */ int zlasr_(char *, char *, char *, integer *, + integer *, doublereal *, doublereal *, doublecomplex *, integer *), zswap_(integer *, doublecomplex *, + integer *, doublecomplex *, integer *), dlaev2_(doublereal *, + doublereal *, doublereal *, doublereal *, doublereal *, + doublereal *, doublereal *); + static integer lendm1, lendp1; + + static integer iscale; + extern /* Subroutine */ int dlascl_(char *, integer *, integer *, + doublereal *, doublereal *, integer *, integer *, doublereal *, + integer *, integer *); + static doublereal safmin; + extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, + doublereal *, doublereal *, doublereal *); + static doublereal safmax; + extern /* Subroutine */ int xerbla_(char *, integer *); + extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *); + extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *, + integer *); + static integer lendsv; + static doublereal ssfmin; + static integer nmaxit, icompz; + static doublereal ssfmax; + extern /* Subroutine */ int zlaset_(char *, integer *, integer *, + doublecomplex *, doublecomplex *, doublecomplex *, integer *); + + +/* + -- LAPACK routine (version 3.0) -- + Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + Courant Institute, Argonne National Lab, and Rice University + September 30, 1994 + + + Purpose + ======= + + ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a + symmetric tridiagonal matrix using the implicit QL or QR method. + The eigenvectors of a full or band complex Hermitian matrix can also + be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this + matrix to tridiagonal form. Arguments ========= - UPLO (input) CHARACTER*1 - = 'U': Upper triangle of A is stored; - = 'L': Lower triangle of A is stored. + COMPZ (input) CHARACTER*1 + = 'N': Compute eigenvalues only. + = 'V': Compute eigenvalues and eigenvectors of the original + Hermitian matrix. On entry, Z must contain the + unitary matrix used to reduce the original matrix + to tridiagonal form. + = 'I': Compute eigenvalues and eigenvectors of the + tridiagonal matrix. Z is initialized to the identity + matrix. N (input) INTEGER - The order of the matrix A. N >= 0. + The order of the matrix. N >= 0. - A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the Hermitian matrix A. If UPLO = 'U', the leading - N-by-N upper triangular part of A contains the upper - triangular part of the matrix A, and the strictly lower - triangular part of A is not referenced. If UPLO = 'L', the - leading N-by-N lower triangular part of A contains the lower - triangular part of the matrix A, and the strictly upper - triangular part of A is not referenced. + D (input/output) DOUBLE PRECISION array, dimension (N) + On entry, the diagonal elements of the tridiagonal matrix. + On exit, if INFO = 0, the eigenvalues in ascending order. - On exit, if INFO = 0, the factor U or L from the Cholesky - factorization A = U**H*U or A = L*L**H. + E (input/output) DOUBLE PRECISION array, dimension (N-1) + On entry, the (n-1) subdiagonal elements of the tridiagonal + matrix. + On exit, E has been destroyed. - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). + Z (input/output) COMPLEX*16 array, dimension (LDZ, N) + On entry, if COMPZ = 'V', then Z contains the unitary + matrix used in the reduction to tridiagonal form. + On exit, if INFO = 0, then if COMPZ = 'V', Z contains the + orthonormal eigenvectors of the original Hermitian matrix, + and if COMPZ = 'I', Z contains the orthonormal eigenvectors + of the symmetric tridiagonal matrix. + If COMPZ = 'N', then Z is not referenced. + + LDZ (input) INTEGER + The leading dimension of the array Z. LDZ >= 1, and if + eigenvectors are desired, then LDZ >= max(1,N). + + WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) + If COMPZ = 'N', then WORK is not referenced. INFO (output) INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value - > 0: if INFO = i, the leading minor of order i is not - positive definite, and the factorization could not be - completed. + > 0: the algorithm has failed to find all the eigenvalues in + a total of 30*N iterations; if INFO = i, then i + elements of E have not converged to zero; on exit, D + and E contain the elements of a symmetric tridiagonal + matrix which is unitarily similar to the original + matrix. ===================================================================== @@ -21158,23 +21502,35 @@ L40: */ /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; + --d__; + --e; + z_dim1 = *ldz; + z_offset = 1 + z_dim1; + z__ -= z_offset; + --work; /* Function Body */ *info = 0; - upper = lsame_(uplo, "U"); - if (! upper && ! lsame_(uplo, "L")) { + + if (lsame_(compz, "N")) { + icompz = 0; + } else if (lsame_(compz, "V")) { + icompz = 1; + } else if (lsame_(compz, "I")) { + icompz = 2; + } else { + icompz = -1; + } + if (icompz < 0) { *info = -1; } else if (*n < 0) { *info = -2; - } else if (*lda < max(1,*n)) { - *info = -4; + } else if (*ldz < 1 || icompz > 0 && *ldz < max(1,*n)) { + *info = -6; } if (*info != 0) { i__1 = -(*info); - xerbla_("ZPOTRF", &i__1); + xerbla_("ZSTEQR", &i__1); return 0; } @@ -21184,2307 +21540,951 @@ L40: return 0; } -/* Determine the block size for this environment. */ - - nb = ilaenv_(&c__1, "ZPOTRF", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( - ftnlen)1); - if ((nb <= 1) || (nb >= *n)) { + if (*n == 1) { + if (icompz == 2) { + i__1 = z_dim1 + 1; + z__[i__1].r = 1., z__[i__1].i = 0.; + } + return 0; + } -/* Use unblocked code. */ +/* Determine the unit roundoff and over/underflow thresholds. */ - zpotf2_(uplo, n, &a[a_offset], lda, info); - } else { - -/* Use blocked code. */ - - if (upper) { - -/* Compute the Cholesky factorization A = U'*U. */ - - i__1 = *n; - i__2 = nb; - for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { + eps = EPSILON; +/* Computing 2nd power */ + d__1 = eps; + eps2 = d__1 * d__1; + safmin = SAFEMINIMUM; + safmax = 1. / safmin; + ssfmax = sqrt(safmax) / 3.; + ssfmin = sqrt(safmin) / eps2; /* - Update and factorize the current diagonal block and test - for non-positive-definiteness. - - Computing MIN + Compute the eigenvalues and eigenvectors of the tridiagonal + matrix. */ - i__3 = nb, i__4 = *n - j + 1; - jb = min(i__3,i__4); - i__3 = j - 1; - zherk_("Upper", "Conjugate transpose", &jb, &i__3, &c_b1294, & - a[j * a_dim1 + 1], lda, &c_b1015, &a[j + j * a_dim1], - lda); - zpotf2_("Upper", &jb, &a[j + j * a_dim1], lda, info); - if (*info != 0) { - goto L30; - } - if (j + jb <= *n) { - -/* Compute the current block row. */ - - i__3 = *n - j - jb + 1; - i__4 = j - 1; - z__1.r = -1., z__1.i = -0.; - zgemm_("Conjugate transpose", "No transpose", &jb, &i__3, - &i__4, &z__1, &a[j * a_dim1 + 1], lda, &a[(j + jb) - * a_dim1 + 1], lda, &c_b60, &a[j + (j + jb) * - a_dim1], lda); - i__3 = *n - j - jb + 1; - ztrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", - &jb, &i__3, &c_b60, &a[j + j * a_dim1], lda, &a[ - j + (j + jb) * a_dim1], lda); - } -/* L10: */ - } - - } else { -/* Compute the Cholesky factorization A = L*L'. */ + if (icompz == 2) { + zlaset_("Full", n, n, &c_b59, &c_b60, &z__[z_offset], ldz); + } - i__2 = *n; - i__1 = nb; - for (j = 1; i__1 < 0 ? j >= i__2 : j <= i__2; j += i__1) { + nmaxit = *n * 30; + jtot = 0; /* - Update and factorize the current diagonal block and test - for non-positive-definiteness. - - Computing MIN + Determine where the matrix splits and choose QL or QR iteration + for each block, according to whether top or bottom diagonal + element is smaller. */ - i__3 = nb, i__4 = *n - j + 1; - jb = min(i__3,i__4); - i__3 = j - 1; - zherk_("Lower", "No transpose", &jb, &i__3, &c_b1294, &a[j + - a_dim1], lda, &c_b1015, &a[j + j * a_dim1], lda); - zpotf2_("Lower", &jb, &a[j + j * a_dim1], lda, info); - if (*info != 0) { - goto L30; - } - if (j + jb <= *n) { -/* Compute the current block column. */ + l1 = 1; + nm1 = *n - 1; - i__3 = *n - j - jb + 1; - i__4 = j - 1; - z__1.r = -1., z__1.i = -0.; - zgemm_("No transpose", "Conjugate transpose", &i__3, &jb, - &i__4, &z__1, &a[j + jb + a_dim1], lda, &a[j + - a_dim1], lda, &c_b60, &a[j + jb + j * a_dim1], - lda); - i__3 = *n - j - jb + 1; - ztrsm_("Right", "Lower", "Conjugate transpose", "Non-unit" - , &i__3, &jb, &c_b60, &a[j + j * a_dim1], lda, &a[ - j + jb + j * a_dim1], lda); - } -/* L20: */ +L10: + if (l1 > *n) { + goto L160; + } + if (l1 > 1) { + e[l1 - 1] = 0.; + } + if (l1 <= nm1) { + i__1 = nm1; + for (m = l1; m <= i__1; ++m) { + tst = (d__1 = e[m], abs(d__1)); + if (tst == 0.) { + goto L30; + } + if (tst <= sqrt((d__1 = d__[m], abs(d__1))) * sqrt((d__2 = d__[m + + 1], abs(d__2))) * eps) { + e[m] = 0.; + goto L30; } +/* L20: */ } } - goto L40; + m = *n; L30: - *info = *info + j - 1; - -L40: - return 0; + l = l1; + lsv = l; + lend = m; + lendsv = lend; + l1 = m + 1; + if (lend == l) { + goto L10; + } -/* End of ZPOTRF */ +/* Scale submatrix in rows and columns L to LEND */ -} /* zpotrf_ */ + i__1 = lend - l + 1; + anorm = dlanst_("I", &i__1, &d__[l], &e[l]); + iscale = 0; + if (anorm == 0.) { + goto L10; + } + if (anorm > ssfmax) { + iscale = 1; + i__1 = lend - l + 1; + dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, + info); + i__1 = lend - l; + dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, + info); + } else if (anorm < ssfmin) { + iscale = 2; + i__1 = lend - l + 1; + dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, + info); + i__1 = lend - l; + dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, + info); + } -/* Subroutine */ int zpotri_(char *uplo, integer *n, doublecomplex *a, - integer *lda, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1; +/* Choose between QL and QR iteration */ - /* Local variables */ - extern logical lsame_(char *, char *); - extern /* Subroutine */ int xerbla_(char *, integer *), zlauum_( - char *, integer *, doublecomplex *, integer *, integer *), - ztrtri_(char *, char *, integer *, doublecomplex *, integer *, - integer *); + if ((d__1 = d__[lend], abs(d__1)) < (d__2 = d__[l], abs(d__2))) { + lend = lsv; + l = lendsv; + } + if (lend > l) { /* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - March 31, 1993 - + QL Iteration - Purpose - ======= + Look for small subdiagonal element. +*/ - ZPOTRI computes the inverse of a complex Hermitian positive definite - matrix A using the Cholesky factorization A = U**H*U or A = L*L**H - computed by ZPOTRF. +L40: + if (l != lend) { + lendm1 = lend - 1; + i__1 = lendm1; + for (m = l; m <= i__1; ++m) { +/* Computing 2nd power */ + d__2 = (d__1 = e[m], abs(d__1)); + tst = d__2 * d__2; + if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m + + 1], abs(d__2)) + safmin) { + goto L60; + } +/* L50: */ + } + } - Arguments - ========= + m = lend; - UPLO (input) CHARACTER*1 - = 'U': Upper triangle of A is stored; - = 'L': Lower triangle of A is stored. +L60: + if (m < lend) { + e[m] = 0.; + } + p = d__[l]; + if (m == l) { + goto L80; + } - N (input) INTEGER - The order of the matrix A. N >= 0. +/* + If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + to compute its eigensystem. +*/ - A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the triangular factor U or L from the Cholesky - factorization A = U**H*U or A = L*L**H, as computed by - ZPOTRF. - On exit, the upper or lower triangle of the (Hermitian) - inverse of A, overwriting the input factor U or L. + if (m == l + 1) { + if (icompz > 0) { + dlaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s); + work[l] = c__; + work[*n - 1 + l] = s; + zlasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], & + z__[l * z_dim1 + 1], ldz); + } else { + dlae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2); + } + d__[l] = rt1; + d__[l + 1] = rt2; + e[l] = 0.; + l += 2; + if (l <= lend) { + goto L40; + } + goto L140; + } - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). + if (jtot == nmaxit) { + goto L140; + } + ++jtot; - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - > 0: if INFO = i, the (i,i) element of the factor U or L is - zero, and the inverse could not be computed. +/* Form shift. */ - ===================================================================== + g = (d__[l + 1] - p) / (e[l] * 2.); + r__ = dlapy2_(&g, &c_b1015); + g = d__[m] - p + e[l] / (g + d_sign(&r__, &g)); + s = 1.; + c__ = 1.; + p = 0.; - Test the input parameters. -*/ +/* Inner loop */ - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; - - /* Function Body */ - *info = 0; - if (! lsame_(uplo, "U") && ! lsame_(uplo, "L")) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*lda < max(1,*n)) { - *info = -4; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZPOTRI", &i__1); - return 0; - } - -/* Quick return if possible */ - - if (*n == 0) { - return 0; - } - -/* Invert the triangular Cholesky factor U or L. */ - - ztrtri_(uplo, "Non-unit", n, &a[a_offset], lda, info); - if (*info > 0) { - return 0; - } - -/* Form inv(U)*inv(U)' or inv(L)'*inv(L). */ - - zlauum_(uplo, n, &a[a_offset], lda, info); - - return 0; - -/* End of ZPOTRI */ - -} /* zpotri_ */ - -/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, - doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, - integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1; - - /* Local variables */ - extern logical lsame_(char *, char *); - static logical upper; - extern /* Subroutine */ int ztrsm_(char *, char *, char *, char *, - integer *, integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *), - xerbla_(char *, integer *); - - -/* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 - - - Purpose - ======= - - ZPOTRS solves a system of linear equations A*X = B with a Hermitian - positive definite matrix A using the Cholesky factorization - A = U**H*U or A = L*L**H computed by ZPOTRF. - - Arguments - ========= - - UPLO (input) CHARACTER*1 - = 'U': Upper triangle of A is stored; - = 'L': Lower triangle of A is stored. - - N (input) INTEGER - The order of the matrix A. N >= 0. - - NRHS (input) INTEGER - The number of right hand sides, i.e., the number of columns - of the matrix B. NRHS >= 0. - - A (input) COMPLEX*16 array, dimension (LDA,N) - The triangular factor U or L from the Cholesky factorization - A = U**H*U or A = L*L**H, as computed by ZPOTRF. - - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). - - B (input/output) COMPLEX*16 array, dimension (LDB,NRHS) - On entry, the right hand side matrix B. - On exit, the solution matrix X. - - LDB (input) INTEGER - The leading dimension of the array B. LDB >= max(1,N). - - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - - ===================================================================== - - - Test the input parameters. -*/ - - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; - b_dim1 = *ldb; - b_offset = 1 + b_dim1; - b -= b_offset; - - /* Function Body */ - *info = 0; - upper = lsame_(uplo, "U"); - if (! upper && ! lsame_(uplo, "L")) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*nrhs < 0) { - *info = -3; - } else if (*lda < max(1,*n)) { - *info = -5; - } else if (*ldb < max(1,*n)) { - *info = -7; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZPOTRS", &i__1); - return 0; - } - -/* Quick return if possible */ - - if ((*n == 0) || (*nrhs == 0)) { - return 0; - } - - if (upper) { - -/* - Solve A*X = B where A = U'*U. - - Solve U'*X = B, overwriting B with X. -*/ - - ztrsm_("Left", "Upper", "Conjugate transpose", "Non-unit", n, nrhs, & - c_b60, &a[a_offset], lda, &b[b_offset], ldb); - -/* Solve U*X = B, overwriting B with X. */ - - ztrsm_("Left", "Upper", "No transpose", "Non-unit", n, nrhs, &c_b60, & - a[a_offset], lda, &b[b_offset], ldb); - } else { - -/* - Solve A*X = B where A = L*L'. - - Solve L*X = B, overwriting B with X. -*/ - - ztrsm_("Left", "Lower", "No transpose", "Non-unit", n, nrhs, &c_b60, & - a[a_offset], lda, &b[b_offset], ldb); - -/* Solve L'*X = B, overwriting B with X. */ - - ztrsm_("Left", "Lower", "Conjugate transpose", "Non-unit", n, nrhs, & - c_b60, &a[a_offset], lda, &b[b_offset], ldb); - } - - return 0; - -/* End of ZPOTRS */ - -} /* zpotrs_ */ - -/* Subroutine */ int zstedc_(char *compz, integer *n, doublereal *d__, - doublereal *e, doublecomplex *z__, integer *ldz, doublecomplex *work, - integer *lwork, doublereal *rwork, integer *lrwork, integer *iwork, - integer *liwork, integer *info) -{ - /* System generated locals */ - integer z_dim1, z_offset, i__1, i__2, i__3, i__4; - doublereal d__1, d__2; - - /* Builtin functions */ - double log(doublereal); - integer pow_ii(integer *, integer *); - double sqrt(doublereal); - - /* Local variables */ - static integer i__, j, k, m; - static doublereal p; - static integer ii, ll, end, lgn; - static doublereal eps, tiny; - extern logical lsame_(char *, char *); - static integer lwmin, start; - extern /* Subroutine */ int zswap_(integer *, doublecomplex *, integer *, - doublecomplex *, integer *), zlaed0_(integer *, integer *, - doublereal *, doublereal *, doublecomplex *, integer *, - doublecomplex *, integer *, doublereal *, integer *, integer *); - - extern /* Subroutine */ int dlascl_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, integer *, doublereal *, - integer *, integer *), dstedc_(char *, integer *, - doublereal *, doublereal *, doublereal *, integer *, doublereal *, - integer *, integer *, integer *, integer *), dlaset_( - char *, integer *, integer *, doublereal *, doublereal *, - doublereal *, integer *), xerbla_(char *, integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *); - extern /* Subroutine */ int dsterf_(integer *, doublereal *, doublereal *, - integer *), zlacrm_(integer *, integer *, doublecomplex *, - integer *, doublereal *, integer *, doublecomplex *, integer *, - doublereal *); - static integer liwmin, icompz; - extern /* Subroutine */ int dsteqr_(char *, integer *, doublereal *, - doublereal *, doublereal *, integer *, doublereal *, integer *), zlacpy_(char *, integer *, integer *, doublecomplex *, - integer *, doublecomplex *, integer *); - static doublereal orgnrm; - static integer lrwmin; - static logical lquery; - static integer smlsiz; - extern /* Subroutine */ int zsteqr_(char *, integer *, doublereal *, - doublereal *, doublecomplex *, integer *, doublereal *, integer *); - - -/* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - June 30, 1999 - - - Purpose - ======= - - ZSTEDC computes all eigenvalues and, optionally, eigenvectors of a - symmetric tridiagonal matrix using the divide and conquer method. - The eigenvectors of a full or band complex Hermitian matrix can also - be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this - matrix to tridiagonal form. - - This code makes very mild assumptions about floating point - arithmetic. It will work on machines with a guard digit in - add/subtract, or on those binary machines without guard digits - which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. - It could conceivably fail on hexadecimal or decimal machines - without guard digits, but we know of none. See DLAED3 for details. - - Arguments - ========= - - COMPZ (input) CHARACTER*1 - = 'N': Compute eigenvalues only. - = 'I': Compute eigenvectors of tridiagonal matrix also. - = 'V': Compute eigenvectors of original Hermitian matrix - also. On entry, Z contains the unitary matrix used - to reduce the original matrix to tridiagonal form. - - N (input) INTEGER - The dimension of the symmetric tridiagonal matrix. N >= 0. - - D (input/output) DOUBLE PRECISION array, dimension (N) - On entry, the diagonal elements of the tridiagonal matrix. - On exit, if INFO = 0, the eigenvalues in ascending order. - - E (input/output) DOUBLE PRECISION array, dimension (N-1) - On entry, the subdiagonal elements of the tridiagonal matrix. - On exit, E has been destroyed. - - Z (input/output) COMPLEX*16 array, dimension (LDZ,N) - On entry, if COMPZ = 'V', then Z contains the unitary - matrix used in the reduction to tridiagonal form. - On exit, if INFO = 0, then if COMPZ = 'V', Z contains the - orthonormal eigenvectors of the original Hermitian matrix, - and if COMPZ = 'I', Z contains the orthonormal eigenvectors - of the symmetric tridiagonal matrix. - If COMPZ = 'N', then Z is not referenced. - - LDZ (input) INTEGER - The leading dimension of the array Z. LDZ >= 1. - If eigenvectors are desired, then LDZ >= max(1,N). - - WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) - On exit, if INFO = 0, WORK(1) returns the optimal LWORK. - - LWORK (input) INTEGER - The dimension of the array WORK. - If COMPZ = 'N' or 'I', or N <= 1, LWORK must be at least 1. - If COMPZ = 'V' and N > 1, LWORK must be at least N*N. - - If LWORK = -1, then a workspace query is assumed; the routine - only calculates the optimal size of the WORK array, returns - this value as the first entry of the WORK array, and no error - message related to LWORK is issued by XERBLA. - - RWORK (workspace/output) DOUBLE PRECISION array, - dimension (LRWORK) - On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK. - - LRWORK (input) INTEGER - The dimension of the array RWORK. - If COMPZ = 'N' or N <= 1, LRWORK must be at least 1. - If COMPZ = 'V' and N > 1, LRWORK must be at least - 1 + 3*N + 2*N*lg N + 3*N**2 , - where lg( N ) = smallest integer k such - that 2**k >= N. - If COMPZ = 'I' and N > 1, LRWORK must be at least - 1 + 4*N + 2*N**2 . - - If LRWORK = -1, then a workspace query is assumed; the - routine only calculates the optimal size of the RWORK array, - returns this value as the first entry of the RWORK array, and - no error message related to LRWORK is issued by XERBLA. - - IWORK (workspace/output) INTEGER array, dimension (LIWORK) - On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. - - LIWORK (input) INTEGER - The dimension of the array IWORK. - If COMPZ = 'N' or N <= 1, LIWORK must be at least 1. - If COMPZ = 'V' or N > 1, LIWORK must be at least - 6 + 6*N + 5*N*lg N. - If COMPZ = 'I' or N > 1, LIWORK must be at least - 3 + 5*N . - - If LIWORK = -1, then a workspace query is assumed; the - routine only calculates the optimal size of the IWORK array, - returns this value as the first entry of the IWORK array, and - no error message related to LIWORK is issued by XERBLA. - - INFO (output) INTEGER - = 0: successful exit. - < 0: if INFO = -i, the i-th argument had an illegal value. - > 0: The algorithm failed to compute an eigenvalue while - working on the submatrix lying in rows and columns - INFO/(N+1) through mod(INFO,N+1). - - Further Details - =============== - - Based on contributions by - Jeff Rutter, Computer Science Division, University of California - at Berkeley, USA - - ===================================================================== - - - Test the input parameters. -*/ - - /* Parameter adjustments */ - --d__; - --e; - z_dim1 = *ldz; - z_offset = 1 + z_dim1; - z__ -= z_offset; - --work; - --rwork; - --iwork; - - /* Function Body */ - *info = 0; - lquery = ((*lwork == -1) || (*lrwork == -1)) || (*liwork == -1); - - if (lsame_(compz, "N")) { - icompz = 0; - } else if (lsame_(compz, "V")) { - icompz = 1; - } else if (lsame_(compz, "I")) { - icompz = 2; - } else { - icompz = -1; - } - if ((*n <= 1) || (icompz <= 0)) { - lwmin = 1; - liwmin = 1; - lrwmin = 1; - } else { - lgn = (integer) (log((doublereal) (*n)) / log(2.)); - if (pow_ii(&c__2, &lgn) < *n) { - ++lgn; - } - if (pow_ii(&c__2, &lgn) < *n) { - ++lgn; - } - if (icompz == 1) { - lwmin = *n * *n; -/* Computing 2nd power */ - i__1 = *n; - lrwmin = *n * 3 + 1 + ((*n) << (1)) * lgn + i__1 * i__1 * 3; - liwmin = *n * 6 + 6 + *n * 5 * lgn; - } else if (icompz == 2) { - lwmin = 1; -/* Computing 2nd power */ - i__1 = *n; - lrwmin = ((*n) << (2)) + 1 + ((i__1 * i__1) << (1)); - liwmin = *n * 5 + 3; - } - } - if (icompz < 0) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if ((*ldz < 1) || (icompz > 0 && *ldz < max(1,*n))) { - *info = -6; - } else if (*lwork < lwmin && ! lquery) { - *info = -8; - } else if (*lrwork < lrwmin && ! lquery) { - *info = -10; - } else if (*liwork < liwmin && ! lquery) { - *info = -12; - } - - if (*info == 0) { - work[1].r = (doublereal) lwmin, work[1].i = 0.; - rwork[1] = (doublereal) lrwmin; - iwork[1] = liwmin; - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZSTEDC", &i__1); - return 0; - } else if (lquery) { - return 0; - } - -/* Quick return if possible */ - - if (*n == 0) { - return 0; - } - if (*n == 1) { - if (icompz != 0) { - i__1 = z_dim1 + 1; - z__[i__1].r = 1., z__[i__1].i = 0.; - } - return 0; - } - - smlsiz = ilaenv_(&c__9, "ZSTEDC", " ", &c__0, &c__0, &c__0, &c__0, ( - ftnlen)6, (ftnlen)1); - -/* - If the following conditional clause is removed, then the routine - will use the Divide and Conquer routine to compute only the - eigenvalues, which requires (3N + 3N**2) real workspace and - (2 + 5N + 2N lg(N)) integer workspace. - Since on many architectures DSTERF is much faster than any other - algorithm for finding eigenvalues only, it is used here - as the default. - - If COMPZ = 'N', use DSTERF to compute the eigenvalues. -*/ - - if (icompz == 0) { - dsterf_(n, &d__[1], &e[1], info); - return 0; - } - -/* - If N is smaller than the minimum divide size (SMLSIZ+1), then - solve the problem with another solver. -*/ - - if (*n <= smlsiz) { - if (icompz == 0) { - dsterf_(n, &d__[1], &e[1], info); - return 0; - } else if (icompz == 2) { - zsteqr_("I", n, &d__[1], &e[1], &z__[z_offset], ldz, &rwork[1], - info); - return 0; - } else { - zsteqr_("V", n, &d__[1], &e[1], &z__[z_offset], ldz, &rwork[1], - info); - return 0; - } - } - -/* If COMPZ = 'I', we simply call DSTEDC instead. */ - - if (icompz == 2) { - dlaset_("Full", n, n, &c_b324, &c_b1015, &rwork[1], n); - ll = *n * *n + 1; - i__1 = *lrwork - ll + 1; - dstedc_("I", n, &d__[1], &e[1], &rwork[1], n, &rwork[ll], &i__1, & - iwork[1], liwork, info); - i__1 = *n; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - i__3 = i__ + j * z_dim1; - i__4 = (j - 1) * *n + i__; - z__[i__3].r = rwork[i__4], z__[i__3].i = 0.; -/* L10: */ - } -/* L20: */ - } - return 0; - } - -/* - From now on, only option left to be handled is COMPZ = 'V', - i.e. ICOMPZ = 1. - - Scale. -*/ - - orgnrm = dlanst_("M", n, &d__[1], &e[1]); - if (orgnrm == 0.) { - return 0; - } - - eps = EPSILON; - - start = 1; - -/* while ( START <= N ) */ - -L30: - if (start <= *n) { - -/* - Let END be the position of the next subdiagonal entry such that - E( END ) <= TINY or END = N if no such subdiagonal exists. The - matrix identified by the elements between START and END - constitutes an independent sub-problem. -*/ - - end = start; -L40: - if (end < *n) { - tiny = eps * sqrt((d__1 = d__[end], abs(d__1))) * sqrt((d__2 = - d__[end + 1], abs(d__2))); - if ((d__1 = e[end], abs(d__1)) > tiny) { - ++end; - goto L40; - } - } - -/* (Sub) Problem determined. Compute its size and solve it. */ - - m = end - start + 1; - if (m > smlsiz) { - *info = smlsiz; - -/* Scale. */ - - orgnrm = dlanst_("M", &m, &d__[start], &e[start]); - dlascl_("G", &c__0, &c__0, &orgnrm, &c_b1015, &m, &c__1, &d__[ - start], &m, info); - i__1 = m - 1; - i__2 = m - 1; - dlascl_("G", &c__0, &c__0, &orgnrm, &c_b1015, &i__1, &c__1, &e[ - start], &i__2, info); - - zlaed0_(n, &m, &d__[start], &e[start], &z__[start * z_dim1 + 1], - ldz, &work[1], n, &rwork[1], &iwork[1], info); - if (*info > 0) { - *info = (*info / (m + 1) + start - 1) * (*n + 1) + *info % (m - + 1) + start - 1; - return 0; - } - -/* Scale back. */ - - dlascl_("G", &c__0, &c__0, &c_b1015, &orgnrm, &m, &c__1, &d__[ - start], &m, info); - - } else { - dsteqr_("I", &m, &d__[start], &e[start], &rwork[1], &m, &rwork[m * - m + 1], info); - zlacrm_(n, &m, &z__[start * z_dim1 + 1], ldz, &rwork[1], &m, & - work[1], n, &rwork[m * m + 1]); - zlacpy_("A", n, &m, &work[1], n, &z__[start * z_dim1 + 1], ldz); - if (*info > 0) { - *info = start * (*n + 1) + end; - return 0; - } - } - - start = end + 1; - goto L30; - } - -/* - endwhile - - If the problem split any number of times, then the eigenvalues - will not be properly ordered. Here we permute the eigenvalues - (and the associated eigenvectors) into ascending order. -*/ - - if (m != *n) { - -/* Use Selection Sort to minimize swaps of eigenvectors */ - - i__1 = *n; - for (ii = 2; ii <= i__1; ++ii) { - i__ = ii - 1; - k = i__; - p = d__[i__]; - i__2 = *n; - for (j = ii; j <= i__2; ++j) { - if (d__[j] < p) { - k = j; - p = d__[j]; - } -/* L50: */ - } - if (k != i__) { - d__[k] = d__[i__]; - d__[i__] = p; - zswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1], - &c__1); - } -/* L60: */ - } - } - - work[1].r = (doublereal) lwmin, work[1].i = 0.; - rwork[1] = (doublereal) lrwmin; - iwork[1] = liwmin; - - return 0; - -/* End of ZSTEDC */ - -} /* zstedc_ */ - -/* Subroutine */ int zsteqr_(char *compz, integer *n, doublereal *d__, - doublereal *e, doublecomplex *z__, integer *ldz, doublereal *work, - integer *info) -{ - /* System generated locals */ - integer z_dim1, z_offset, i__1, i__2; - doublereal d__1, d__2; - - /* Builtin functions */ - double sqrt(doublereal), d_sign(doublereal *, doublereal *); - - /* Local variables */ - static doublereal b, c__, f, g; - static integer i__, j, k, l, m; - static doublereal p, r__, s; - static integer l1, ii, mm, lm1, mm1, nm1; - static doublereal rt1, rt2, eps; - static integer lsv; - static doublereal tst, eps2; - static integer lend, jtot; - extern /* Subroutine */ int dlae2_(doublereal *, doublereal *, doublereal - *, doublereal *, doublereal *); - extern logical lsame_(char *, char *); - static doublereal anorm; - extern /* Subroutine */ int zlasr_(char *, char *, char *, integer *, - integer *, doublereal *, doublereal *, doublecomplex *, integer *), zswap_(integer *, doublecomplex *, - integer *, doublecomplex *, integer *), dlaev2_(doublereal *, - doublereal *, doublereal *, doublereal *, doublereal *, - doublereal *, doublereal *); - static integer lendm1, lendp1; - - static integer iscale; - extern /* Subroutine */ int dlascl_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, integer *, doublereal *, - integer *, integer *); - static doublereal safmin; - extern /* Subroutine */ int dlartg_(doublereal *, doublereal *, - doublereal *, doublereal *, doublereal *); - static doublereal safmax; - extern /* Subroutine */ int xerbla_(char *, integer *); - extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *); - extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *, - integer *); - static integer lendsv; - static doublereal ssfmin; - static integer nmaxit, icompz; - static doublereal ssfmax; - extern /* Subroutine */ int zlaset_(char *, integer *, integer *, - doublecomplex *, doublecomplex *, doublecomplex *, integer *); - - -/* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 - - - Purpose - ======= - - ZSTEQR computes all eigenvalues and, optionally, eigenvectors of a - symmetric tridiagonal matrix using the implicit QL or QR method. - The eigenvectors of a full or band complex Hermitian matrix can also - be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this - matrix to tridiagonal form. - - Arguments - ========= - - COMPZ (input) CHARACTER*1 - = 'N': Compute eigenvalues only. - = 'V': Compute eigenvalues and eigenvectors of the original - Hermitian matrix. On entry, Z must contain the - unitary matrix used to reduce the original matrix - to tridiagonal form. - = 'I': Compute eigenvalues and eigenvectors of the - tridiagonal matrix. Z is initialized to the identity - matrix. - - N (input) INTEGER - The order of the matrix. N >= 0. - - D (input/output) DOUBLE PRECISION array, dimension (N) - On entry, the diagonal elements of the tridiagonal matrix. - On exit, if INFO = 0, the eigenvalues in ascending order. - - E (input/output) DOUBLE PRECISION array, dimension (N-1) - On entry, the (n-1) subdiagonal elements of the tridiagonal - matrix. - On exit, E has been destroyed. - - Z (input/output) COMPLEX*16 array, dimension (LDZ, N) - On entry, if COMPZ = 'V', then Z contains the unitary - matrix used in the reduction to tridiagonal form. - On exit, if INFO = 0, then if COMPZ = 'V', Z contains the - orthonormal eigenvectors of the original Hermitian matrix, - and if COMPZ = 'I', Z contains the orthonormal eigenvectors - of the symmetric tridiagonal matrix. - If COMPZ = 'N', then Z is not referenced. - - LDZ (input) INTEGER - The leading dimension of the array Z. LDZ >= 1, and if - eigenvectors are desired, then LDZ >= max(1,N). - - WORK (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2)) - If COMPZ = 'N', then WORK is not referenced. - - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - > 0: the algorithm has failed to find all the eigenvalues in - a total of 30*N iterations; if INFO = i, then i - elements of E have not converged to zero; on exit, D - and E contain the elements of a symmetric tridiagonal - matrix which is unitarily similar to the original - matrix. - - ===================================================================== - - - Test the input parameters. -*/ - - /* Parameter adjustments */ - --d__; - --e; - z_dim1 = *ldz; - z_offset = 1 + z_dim1; - z__ -= z_offset; - --work; - - /* Function Body */ - *info = 0; - - if (lsame_(compz, "N")) { - icompz = 0; - } else if (lsame_(compz, "V")) { - icompz = 1; - } else if (lsame_(compz, "I")) { - icompz = 2; - } else { - icompz = -1; - } - if (icompz < 0) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if ((*ldz < 1) || (icompz > 0 && *ldz < max(1,*n))) { - *info = -6; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZSTEQR", &i__1); - return 0; - } - -/* Quick return if possible */ - - if (*n == 0) { - return 0; - } - - if (*n == 1) { - if (icompz == 2) { - i__1 = z_dim1 + 1; - z__[i__1].r = 1., z__[i__1].i = 0.; - } - return 0; - } - -/* Determine the unit roundoff and over/underflow thresholds. */ - - eps = EPSILON; -/* Computing 2nd power */ - d__1 = eps; - eps2 = d__1 * d__1; - safmin = SAFEMINIMUM; - safmax = 1. / safmin; - ssfmax = sqrt(safmax) / 3.; - ssfmin = sqrt(safmin) / eps2; - -/* - Compute the eigenvalues and eigenvectors of the tridiagonal - matrix. -*/ - - if (icompz == 2) { - zlaset_("Full", n, n, &c_b59, &c_b60, &z__[z_offset], ldz); - } - - nmaxit = *n * 30; - jtot = 0; - -/* - Determine where the matrix splits and choose QL or QR iteration - for each block, according to whether top or bottom diagonal - element is smaller. -*/ - - l1 = 1; - nm1 = *n - 1; - -L10: - if (l1 > *n) { - goto L160; - } - if (l1 > 1) { - e[l1 - 1] = 0.; - } - if (l1 <= nm1) { - i__1 = nm1; - for (m = l1; m <= i__1; ++m) { - tst = (d__1 = e[m], abs(d__1)); - if (tst == 0.) { - goto L30; - } - if (tst <= sqrt((d__1 = d__[m], abs(d__1))) * sqrt((d__2 = d__[m - + 1], abs(d__2))) * eps) { - e[m] = 0.; - goto L30; - } -/* L20: */ - } - } - m = *n; - -L30: - l = l1; - lsv = l; - lend = m; - lendsv = lend; - l1 = m + 1; - if (lend == l) { - goto L10; - } - -/* Scale submatrix in rows and columns L to LEND */ - - i__1 = lend - l + 1; - anorm = dlanst_("I", &i__1, &d__[l], &e[l]); - iscale = 0; - if (anorm == 0.) { - goto L10; - } - if (anorm > ssfmax) { - iscale = 1; - i__1 = lend - l + 1; - dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n, - info); - i__1 = lend - l; - dlascl_("G", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n, - info); - } else if (anorm < ssfmin) { - iscale = 2; - i__1 = lend - l + 1; - dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n, - info); - i__1 = lend - l; - dlascl_("G", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n, - info); - } - -/* Choose between QL and QR iteration */ - - if ((d__1 = d__[lend], abs(d__1)) < (d__2 = d__[l], abs(d__2))) { - lend = lsv; - l = lendsv; - } - - if (lend > l) { - -/* - QL Iteration - - Look for small subdiagonal element. -*/ - -L40: - if (l != lend) { - lendm1 = lend - 1; - i__1 = lendm1; - for (m = l; m <= i__1; ++m) { -/* Computing 2nd power */ - d__2 = (d__1 = e[m], abs(d__1)); - tst = d__2 * d__2; - if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m - + 1], abs(d__2)) + safmin) { - goto L60; - } -/* L50: */ - } - } - - m = lend; - -L60: - if (m < lend) { - e[m] = 0.; - } - p = d__[l]; - if (m == l) { - goto L80; - } - -/* - If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 - to compute its eigensystem. -*/ - - if (m == l + 1) { - if (icompz > 0) { - dlaev2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s); - work[l] = c__; - work[*n - 1 + l] = s; - zlasr_("R", "V", "B", n, &c__2, &work[l], &work[*n - 1 + l], & - z__[l * z_dim1 + 1], ldz); - } else { - dlae2_(&d__[l], &e[l], &d__[l + 1], &rt1, &rt2); - } - d__[l] = rt1; - d__[l + 1] = rt2; - e[l] = 0.; - l += 2; - if (l <= lend) { - goto L40; - } - goto L140; - } - - if (jtot == nmaxit) { - goto L140; - } - ++jtot; - -/* Form shift. */ - - g = (d__[l + 1] - p) / (e[l] * 2.); - r__ = dlapy2_(&g, &c_b1015); - g = d__[m] - p + e[l] / (g + d_sign(&r__, &g)); - - s = 1.; - c__ = 1.; - p = 0.; - -/* Inner loop */ - - mm1 = m - 1; - i__1 = l; - for (i__ = mm1; i__ >= i__1; --i__) { - f = s * e[i__]; - b = c__ * e[i__]; - dlartg_(&g, &f, &c__, &s, &r__); - if (i__ != m - 1) { - e[i__ + 1] = r__; - } - g = d__[i__ + 1] - p; - r__ = (d__[i__] - g) * s + c__ * 2. * b; - p = s * r__; - d__[i__ + 1] = g + p; - g = c__ * r__ - b; - -/* If eigenvectors are desired, then save rotations. */ - - if (icompz > 0) { - work[i__] = c__; - work[*n - 1 + i__] = -s; - } - -/* L70: */ - } - -/* If eigenvectors are desired, then apply saved rotations. */ - - if (icompz > 0) { - mm = m - l + 1; - zlasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l - * z_dim1 + 1], ldz); - } - - d__[l] -= p; - e[l] = g; - goto L40; - -/* Eigenvalue found. */ - -L80: - d__[l] = p; - - ++l; - if (l <= lend) { - goto L40; - } - goto L140; - - } else { - -/* - QR Iteration - - Look for small superdiagonal element. -*/ - -L90: - if (l != lend) { - lendp1 = lend + 1; - i__1 = lendp1; - for (m = l; m >= i__1; --m) { -/* Computing 2nd power */ - d__2 = (d__1 = e[m - 1], abs(d__1)); - tst = d__2 * d__2; - if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m - - 1], abs(d__2)) + safmin) { - goto L110; - } -/* L100: */ - } - } - - m = lend; - -L110: - if (m > lend) { - e[m - 1] = 0.; - } - p = d__[l]; - if (m == l) { - goto L130; - } - -/* - If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 - to compute its eigensystem. -*/ - - if (m == l - 1) { - if (icompz > 0) { - dlaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s) - ; - work[m] = c__; - work[*n - 1 + m] = s; - zlasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], & - z__[(l - 1) * z_dim1 + 1], ldz); - } else { - dlae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2); - } - d__[l - 1] = rt1; - d__[l] = rt2; - e[l - 1] = 0.; - l += -2; - if (l >= lend) { - goto L90; - } - goto L140; - } - - if (jtot == nmaxit) { - goto L140; - } - ++jtot; - -/* Form shift. */ - - g = (d__[l - 1] - p) / (e[l - 1] * 2.); - r__ = dlapy2_(&g, &c_b1015); - g = d__[m] - p + e[l - 1] / (g + d_sign(&r__, &g)); - - s = 1.; - c__ = 1.; - p = 0.; - -/* Inner loop */ - - lm1 = l - 1; - i__1 = lm1; - for (i__ = m; i__ <= i__1; ++i__) { - f = s * e[i__]; - b = c__ * e[i__]; - dlartg_(&g, &f, &c__, &s, &r__); - if (i__ != m) { - e[i__ - 1] = r__; - } - g = d__[i__] - p; - r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b; - p = s * r__; - d__[i__] = g + p; - g = c__ * r__ - b; - -/* If eigenvectors are desired, then save rotations. */ - - if (icompz > 0) { - work[i__] = c__; - work[*n - 1 + i__] = s; - } - -/* L120: */ - } - -/* If eigenvectors are desired, then apply saved rotations. */ - - if (icompz > 0) { - mm = l - m + 1; - zlasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m - * z_dim1 + 1], ldz); - } - - d__[l] -= p; - e[lm1] = g; - goto L90; - -/* Eigenvalue found. */ - -L130: - d__[l] = p; - - --l; - if (l >= lend) { - goto L90; - } - goto L140; - - } - -/* Undo scaling if necessary */ - -L140: - if (iscale == 1) { - i__1 = lendsv - lsv + 1; - dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], - n, info); - i__1 = lendsv - lsv; - dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n, - info); - } else if (iscale == 2) { - i__1 = lendsv - lsv + 1; - dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], - n, info); - i__1 = lendsv - lsv; - dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n, - info); - } - -/* - Check for no convergence to an eigenvalue after a total - of N*MAXIT iterations. -*/ - - if (jtot == nmaxit) { - i__1 = *n - 1; - for (i__ = 1; i__ <= i__1; ++i__) { - if (e[i__] != 0.) { - ++(*info); - } -/* L150: */ - } - return 0; - } - goto L10; - -/* Order eigenvalues and eigenvectors. */ - -L160: - if (icompz == 0) { - -/* Use Quick Sort */ - - dlasrt_("I", n, &d__[1], info); - - } else { - -/* Use Selection Sort to minimize swaps of eigenvectors */ - - i__1 = *n; - for (ii = 2; ii <= i__1; ++ii) { - i__ = ii - 1; - k = i__; - p = d__[i__]; - i__2 = *n; - for (j = ii; j <= i__2; ++j) { - if (d__[j] < p) { - k = j; - p = d__[j]; - } -/* L170: */ - } - if (k != i__) { - d__[k] = d__[i__]; - d__[i__] = p; - zswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1], - &c__1); - } -/* L180: */ - } - } - return 0; - -/* End of ZSTEQR */ - -} /* zsteqr_ */ - -/* Subroutine */ int ztrevc_(char *side, char *howmny, logical *select, - integer *n, doublecomplex *t, integer *ldt, doublecomplex *vl, - integer *ldvl, doublecomplex *vr, integer *ldvr, integer *mm, integer - *m, doublecomplex *work, doublereal *rwork, integer *info) -{ - /* System generated locals */ - integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, - i__2, i__3, i__4, i__5; - doublereal d__1, d__2, d__3; - doublecomplex z__1, z__2; - - /* Builtin functions */ - double d_imag(doublecomplex *); - void d_cnjg(doublecomplex *, doublecomplex *); - - /* Local variables */ - static integer i__, j, k, ii, ki, is; - static doublereal ulp; - static logical allv; - static doublereal unfl, ovfl, smin; - static logical over; - static doublereal scale; - extern logical lsame_(char *, char *); - static doublereal remax; - static logical leftv, bothv; - extern /* Subroutine */ int zgemv_(char *, integer *, integer *, - doublecomplex *, doublecomplex *, integer *, doublecomplex *, - integer *, doublecomplex *, doublecomplex *, integer *); - static logical somev; - extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, - doublecomplex *, integer *), dlabad_(doublereal *, doublereal *); - - extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_( - integer *, doublereal *, doublecomplex *, integer *); - extern integer izamax_(integer *, doublecomplex *, integer *); - static logical rightv; - extern doublereal dzasum_(integer *, doublecomplex *, integer *); - static doublereal smlnum; - extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *, - integer *, doublecomplex *, integer *, doublecomplex *, - doublereal *, doublereal *, integer *); - - -/* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - June 30, 1999 - - - Purpose - ======= - - ZTREVC computes some or all of the right and/or left eigenvectors of - a complex upper triangular matrix T. - - The right eigenvector x and the left eigenvector y of T corresponding - to an eigenvalue w are defined by: - - T*x = w*x, y'*T = w*y' - - where y' denotes the conjugate transpose of the vector y. - - If all eigenvectors are requested, the routine may either return the - matrices X and/or Y of right or left eigenvectors of T, or the - products Q*X and/or Q*Y, where Q is an input unitary - matrix. If T was obtained from the Schur factorization of an - original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of - right or left eigenvectors of A. - - Arguments - ========= - - SIDE (input) CHARACTER*1 - = 'R': compute right eigenvectors only; - = 'L': compute left eigenvectors only; - = 'B': compute both right and left eigenvectors. - - HOWMNY (input) CHARACTER*1 - = 'A': compute all right and/or left eigenvectors; - = 'B': compute all right and/or left eigenvectors, - and backtransform them using the input matrices - supplied in VR and/or VL; - = 'S': compute selected right and/or left eigenvectors, - specified by the logical array SELECT. - - SELECT (input) LOGICAL array, dimension (N) - If HOWMNY = 'S', SELECT specifies the eigenvectors to be - computed. - If HOWMNY = 'A' or 'B', SELECT is not referenced. - To select the eigenvector corresponding to the j-th - eigenvalue, SELECT(j) must be set to .TRUE.. - - N (input) INTEGER - The order of the matrix T. N >= 0. - - T (input/output) COMPLEX*16 array, dimension (LDT,N) - The upper triangular matrix T. T is modified, but restored - on exit. - - LDT (input) INTEGER - The leading dimension of the array T. LDT >= max(1,N). - - VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) - On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must - contain an N-by-N matrix Q (usually the unitary matrix Q of - Schur vectors returned by ZHSEQR). - On exit, if SIDE = 'L' or 'B', VL contains: - if HOWMNY = 'A', the matrix Y of left eigenvectors of T; - VL is lower triangular. The i-th column - VL(i) of VL is the eigenvector corresponding - to T(i,i). - if HOWMNY = 'B', the matrix Q*Y; - if HOWMNY = 'S', the left eigenvectors of T specified by - SELECT, stored consecutively in the columns - of VL, in the same order as their - eigenvalues. - If SIDE = 'R', VL is not referenced. - - LDVL (input) INTEGER - The leading dimension of the array VL. LDVL >= max(1,N) if - SIDE = 'L' or 'B'; LDVL >= 1 otherwise. - - VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) - On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must - contain an N-by-N matrix Q (usually the unitary matrix Q of - Schur vectors returned by ZHSEQR). - On exit, if SIDE = 'R' or 'B', VR contains: - if HOWMNY = 'A', the matrix X of right eigenvectors of T; - VR is upper triangular. The i-th column - VR(i) of VR is the eigenvector corresponding - to T(i,i). - if HOWMNY = 'B', the matrix Q*X; - if HOWMNY = 'S', the right eigenvectors of T specified by - SELECT, stored consecutively in the columns - of VR, in the same order as their - eigenvalues. - If SIDE = 'L', VR is not referenced. - - LDVR (input) INTEGER - The leading dimension of the array VR. LDVR >= max(1,N) if - SIDE = 'R' or 'B'; LDVR >= 1 otherwise. - - MM (input) INTEGER - The number of columns in the arrays VL and/or VR. MM >= M. - - M (output) INTEGER - The number of columns in the arrays VL and/or VR actually - used to store the eigenvectors. If HOWMNY = 'A' or 'B', M - is set to N. Each selected eigenvector occupies one - column. - - WORK (workspace) COMPLEX*16 array, dimension (2*N) - - RWORK (workspace) DOUBLE PRECISION array, dimension (N) - - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - - Further Details - =============== - - The algorithm used in this program is basically backward (forward) - substitution, with scaling to make the code robust against - possible overflow. - - Each eigenvector is normalized so that the element of largest - magnitude has magnitude 1; here the magnitude of a complex number - (x,y) is taken to be |x| + |y|. - - ===================================================================== - - - Decode and test the input parameters -*/ - - /* Parameter adjustments */ - --select; - t_dim1 = *ldt; - t_offset = 1 + t_dim1; - t -= t_offset; - vl_dim1 = *ldvl; - vl_offset = 1 + vl_dim1; - vl -= vl_offset; - vr_dim1 = *ldvr; - vr_offset = 1 + vr_dim1; - vr -= vr_offset; - --work; - --rwork; - - /* Function Body */ - bothv = lsame_(side, "B"); - rightv = (lsame_(side, "R")) || (bothv); - leftv = (lsame_(side, "L")) || (bothv); - - allv = lsame_(howmny, "A"); - over = lsame_(howmny, "B"); - somev = lsame_(howmny, "S"); + mm1 = m - 1; + i__1 = l; + for (i__ = mm1; i__ >= i__1; --i__) { + f = s * e[i__]; + b = c__ * e[i__]; + dlartg_(&g, &f, &c__, &s, &r__); + if (i__ != m - 1) { + e[i__ + 1] = r__; + } + g = d__[i__ + 1] - p; + r__ = (d__[i__] - g) * s + c__ * 2. * b; + p = s * r__; + d__[i__ + 1] = g + p; + g = c__ * r__ - b; -/* - Set M to the number of columns required to store the selected - eigenvectors. -*/ +/* If eigenvectors are desired, then save rotations. */ - if (somev) { - *m = 0; - i__1 = *n; - for (j = 1; j <= i__1; ++j) { - if (select[j]) { - ++(*m); + if (icompz > 0) { + work[i__] = c__; + work[*n - 1 + i__] = -s; } -/* L10: */ + +/* L70: */ } - } else { - *m = *n; - } - *info = 0; - if (! rightv && ! leftv) { - *info = -1; - } else if (! allv && ! over && ! somev) { - *info = -2; - } else if (*n < 0) { - *info = -4; - } else if (*ldt < max(1,*n)) { - *info = -6; - } else if ((*ldvl < 1) || (leftv && *ldvl < *n)) { - *info = -8; - } else if ((*ldvr < 1) || (rightv && *ldvr < *n)) { - *info = -10; - } else if (*mm < *m) { - *info = -11; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZTREVC", &i__1); - return 0; - } +/* If eigenvectors are desired, then apply saved rotations. */ -/* Quick return if possible. */ + if (icompz > 0) { + mm = m - l + 1; + zlasr_("R", "V", "B", n, &mm, &work[l], &work[*n - 1 + l], &z__[l + * z_dim1 + 1], ldz); + } - if (*n == 0) { - return 0; - } + d__[l] -= p; + e[l] = g; + goto L40; -/* Set the constants to control overflow. */ +/* Eigenvalue found. */ - unfl = SAFEMINIMUM; - ovfl = 1. / unfl; - dlabad_(&unfl, &ovfl); - ulp = PRECISION; - smlnum = unfl * (*n / ulp); +L80: + d__[l] = p; -/* Store the diagonal elements of T in working array WORK. */ + ++l; + if (l <= lend) { + goto L40; + } + goto L140; - i__1 = *n; - for (i__ = 1; i__ <= i__1; ++i__) { - i__2 = i__ + *n; - i__3 = i__ + i__ * t_dim1; - work[i__2].r = t[i__3].r, work[i__2].i = t[i__3].i; -/* L20: */ - } + } else { /* - Compute 1-norm of each column of strictly upper triangular - part of T to control overflow in triangular solver. + QR Iteration + + Look for small superdiagonal element. */ - rwork[1] = 0.; - i__1 = *n; - for (j = 2; j <= i__1; ++j) { - i__2 = j - 1; - rwork[j] = dzasum_(&i__2, &t[j * t_dim1 + 1], &c__1); -/* L30: */ - } +L90: + if (l != lend) { + lendp1 = lend + 1; + i__1 = lendp1; + for (m = l; m >= i__1; --m) { +/* Computing 2nd power */ + d__2 = (d__1 = e[m - 1], abs(d__1)); + tst = d__2 * d__2; + if (tst <= eps2 * (d__1 = d__[m], abs(d__1)) * (d__2 = d__[m + - 1], abs(d__2)) + safmin) { + goto L110; + } +/* L100: */ + } + } - if (rightv) { + m = lend; -/* Compute right eigenvectors. */ +L110: + if (m > lend) { + e[m - 1] = 0.; + } + p = d__[l]; + if (m == l) { + goto L130; + } - is = *m; - for (ki = *n; ki >= 1; --ki) { +/* + If remaining matrix is 2-by-2, use DLAE2 or SLAEV2 + to compute its eigensystem. +*/ - if (somev) { - if (! select[ki]) { - goto L80; - } + if (m == l - 1) { + if (icompz > 0) { + dlaev2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s) + ; + work[m] = c__; + work[*n - 1 + m] = s; + zlasr_("R", "V", "F", n, &c__2, &work[m], &work[*n - 1 + m], & + z__[(l - 1) * z_dim1 + 1], ldz); + } else { + dlae2_(&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2); } -/* Computing MAX */ - i__1 = ki + ki * t_dim1; - d__3 = ulp * ((d__1 = t[i__1].r, abs(d__1)) + (d__2 = d_imag(&t[ - ki + ki * t_dim1]), abs(d__2))); - smin = max(d__3,smlnum); + d__[l - 1] = rt1; + d__[l] = rt2; + e[l - 1] = 0.; + l += -2; + if (l >= lend) { + goto L90; + } + goto L140; + } - work[1].r = 1., work[1].i = 0.; + if (jtot == nmaxit) { + goto L140; + } + ++jtot; -/* Form right-hand side. */ +/* Form shift. */ - i__1 = ki - 1; - for (k = 1; k <= i__1; ++k) { - i__2 = k; - i__3 = k + ki * t_dim1; - z__1.r = -t[i__3].r, z__1.i = -t[i__3].i; - work[i__2].r = z__1.r, work[i__2].i = z__1.i; -/* L40: */ - } + g = (d__[l - 1] - p) / (e[l - 1] * 2.); + r__ = dlapy2_(&g, &c_b1015); + g = d__[m] - p + e[l - 1] / (g + d_sign(&r__, &g)); -/* - Solve the triangular system: - (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. -*/ + s = 1.; + c__ = 1.; + p = 0.; - i__1 = ki - 1; - for (k = 1; k <= i__1; ++k) { - i__2 = k + k * t_dim1; - i__3 = k + k * t_dim1; - i__4 = ki + ki * t_dim1; - z__1.r = t[i__3].r - t[i__4].r, z__1.i = t[i__3].i - t[i__4] - .i; - t[i__2].r = z__1.r, t[i__2].i = z__1.i; - i__2 = k + k * t_dim1; - if ((d__1 = t[i__2].r, abs(d__1)) + (d__2 = d_imag(&t[k + k * - t_dim1]), abs(d__2)) < smin) { - i__3 = k + k * t_dim1; - t[i__3].r = smin, t[i__3].i = 0.; - } -/* L50: */ - } +/* Inner loop */ - if (ki > 1) { - i__1 = ki - 1; - zlatrs_("Upper", "No transpose", "Non-unit", "Y", &i__1, &t[ - t_offset], ldt, &work[1], &scale, &rwork[1], info); - i__1 = ki; - work[i__1].r = scale, work[i__1].i = 0.; + lm1 = l - 1; + i__1 = lm1; + for (i__ = m; i__ <= i__1; ++i__) { + f = s * e[i__]; + b = c__ * e[i__]; + dlartg_(&g, &f, &c__, &s, &r__); + if (i__ != m) { + e[i__ - 1] = r__; } + g = d__[i__] - p; + r__ = (d__[i__ + 1] - g) * s + c__ * 2. * b; + p = s * r__; + d__[i__] = g + p; + g = c__ * r__ - b; -/* Copy the vector x or Q*x to VR and normalize. */ +/* If eigenvectors are desired, then save rotations. */ - if (! over) { - zcopy_(&ki, &work[1], &c__1, &vr[is * vr_dim1 + 1], &c__1); + if (icompz > 0) { + work[i__] = c__; + work[*n - 1 + i__] = s; + } - ii = izamax_(&ki, &vr[is * vr_dim1 + 1], &c__1); - i__1 = ii + is * vr_dim1; - remax = 1. / ((d__1 = vr[i__1].r, abs(d__1)) + (d__2 = d_imag( - &vr[ii + is * vr_dim1]), abs(d__2))); - zdscal_(&ki, &remax, &vr[is * vr_dim1 + 1], &c__1); +/* L120: */ + } - i__1 = *n; - for (k = ki + 1; k <= i__1; ++k) { - i__2 = k + is * vr_dim1; - vr[i__2].r = 0., vr[i__2].i = 0.; -/* L60: */ - } - } else { - if (ki > 1) { - i__1 = ki - 1; - z__1.r = scale, z__1.i = 0.; - zgemv_("N", n, &i__1, &c_b60, &vr[vr_offset], ldvr, &work[ - 1], &c__1, &z__1, &vr[ki * vr_dim1 + 1], &c__1); - } +/* If eigenvectors are desired, then apply saved rotations. */ + + if (icompz > 0) { + mm = l - m + 1; + zlasr_("R", "V", "F", n, &mm, &work[m], &work[*n - 1 + m], &z__[m + * z_dim1 + 1], ldz); + } - ii = izamax_(n, &vr[ki * vr_dim1 + 1], &c__1); - i__1 = ii + ki * vr_dim1; - remax = 1. / ((d__1 = vr[i__1].r, abs(d__1)) + (d__2 = d_imag( - &vr[ii + ki * vr_dim1]), abs(d__2))); - zdscal_(n, &remax, &vr[ki * vr_dim1 + 1], &c__1); - } + d__[l] -= p; + e[lm1] = g; + goto L90; -/* Set back the original diagonal elements of T. */ +/* Eigenvalue found. */ - i__1 = ki - 1; - for (k = 1; k <= i__1; ++k) { - i__2 = k + k * t_dim1; - i__3 = k + *n; - t[i__2].r = work[i__3].r, t[i__2].i = work[i__3].i; -/* L70: */ - } +L130: + d__[l] = p; - --is; -L80: - ; + --l; + if (l >= lend) { + goto L90; } + goto L140; + } - if (leftv) { +/* Undo scaling if necessary */ -/* Compute left eigenvectors. */ +L140: + if (iscale == 1) { + i__1 = lendsv - lsv + 1; + dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv], + n, info); + i__1 = lendsv - lsv; + dlascl_("G", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n, + info); + } else if (iscale == 2) { + i__1 = lendsv - lsv + 1; + dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv], + n, info); + i__1 = lendsv - lsv; + dlascl_("G", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n, + info); + } - is = 1; - i__1 = *n; - for (ki = 1; ki <= i__1; ++ki) { +/* + Check for no convergence to an eigenvalue after a total + of N*MAXIT iterations. +*/ - if (somev) { - if (! select[ki]) { - goto L130; - } + if (jtot == nmaxit) { + i__1 = *n - 1; + for (i__ = 1; i__ <= i__1; ++i__) { + if (e[i__] != 0.) { + ++(*info); } -/* Computing MAX */ - i__2 = ki + ki * t_dim1; - d__3 = ulp * ((d__1 = t[i__2].r, abs(d__1)) + (d__2 = d_imag(&t[ - ki + ki * t_dim1]), abs(d__2))); - smin = max(d__3,smlnum); +/* L150: */ + } + return 0; + } + goto L10; - i__2 = *n; - work[i__2].r = 1., work[i__2].i = 0.; +/* Order eigenvalues and eigenvectors. */ -/* Form right-hand side. */ +L160: + if (icompz == 0) { - i__2 = *n; - for (k = ki + 1; k <= i__2; ++k) { - i__3 = k; - d_cnjg(&z__2, &t[ki + k * t_dim1]); - z__1.r = -z__2.r, z__1.i = -z__2.i; - work[i__3].r = z__1.r, work[i__3].i = z__1.i; -/* L90: */ - } +/* Use Quick Sort */ -/* - Solve the triangular system: - (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. -*/ + dlasrt_("I", n, &d__[1], info); + + } else { + +/* Use Selection Sort to minimize swaps of eigenvectors */ + i__1 = *n; + for (ii = 2; ii <= i__1; ++ii) { + i__ = ii - 1; + k = i__; + p = d__[i__]; i__2 = *n; - for (k = ki + 1; k <= i__2; ++k) { - i__3 = k + k * t_dim1; - i__4 = k + k * t_dim1; - i__5 = ki + ki * t_dim1; - z__1.r = t[i__4].r - t[i__5].r, z__1.i = t[i__4].i - t[i__5] - .i; - t[i__3].r = z__1.r, t[i__3].i = z__1.i; - i__3 = k + k * t_dim1; - if ((d__1 = t[i__3].r, abs(d__1)) + (d__2 = d_imag(&t[k + k * - t_dim1]), abs(d__2)) < smin) { - i__4 = k + k * t_dim1; - t[i__4].r = smin, t[i__4].i = 0.; + for (j = ii; j <= i__2; ++j) { + if (d__[j] < p) { + k = j; + p = d__[j]; } -/* L100: */ +/* L170: */ } - - if (ki < *n) { - i__2 = *n - ki; - zlatrs_("Upper", "Conjugate transpose", "Non-unit", "Y", & - i__2, &t[ki + 1 + (ki + 1) * t_dim1], ldt, &work[ki + - 1], &scale, &rwork[1], info); - i__2 = ki; - work[i__2].r = scale, work[i__2].i = 0.; + if (k != i__) { + d__[k] = d__[i__]; + d__[i__] = p; + zswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[k * z_dim1 + 1], + &c__1); } +/* L180: */ + } + } + return 0; -/* Copy the vector x or Q*x to VL and normalize. */ +/* End of ZSTEQR */ - if (! over) { - i__2 = *n - ki + 1; - zcopy_(&i__2, &work[ki], &c__1, &vl[ki + is * vl_dim1], &c__1) - ; +} /* zsteqr_ */ - i__2 = *n - ki + 1; - ii = izamax_(&i__2, &vl[ki + is * vl_dim1], &c__1) + ki - 1; - i__2 = ii + is * vl_dim1; - remax = 1. / ((d__1 = vl[i__2].r, abs(d__1)) + (d__2 = d_imag( - &vl[ii + is * vl_dim1]), abs(d__2))); - i__2 = *n - ki + 1; - zdscal_(&i__2, &remax, &vl[ki + is * vl_dim1], &c__1); +/* Subroutine */ int ztrevc_(char *side, char *howmny, logical *select, + integer *n, doublecomplex *t, integer *ldt, doublecomplex *vl, + integer *ldvl, doublecomplex *vr, integer *ldvr, integer *mm, integer + *m, doublecomplex *work, doublereal *rwork, integer *info) +{ + /* System generated locals */ + integer t_dim1, t_offset, vl_dim1, vl_offset, vr_dim1, vr_offset, i__1, + i__2, i__3, i__4, i__5; + doublereal d__1, d__2, d__3; + doublecomplex z__1, z__2; - i__2 = ki - 1; - for (k = 1; k <= i__2; ++k) { - i__3 = k + is * vl_dim1; - vl[i__3].r = 0., vl[i__3].i = 0.; -/* L110: */ - } - } else { - if (ki < *n) { - i__2 = *n - ki; - z__1.r = scale, z__1.i = 0.; - zgemv_("N", n, &i__2, &c_b60, &vl[(ki + 1) * vl_dim1 + 1], - ldvl, &work[ki + 1], &c__1, &z__1, &vl[ki * - vl_dim1 + 1], &c__1); - } + /* Builtin functions */ + double d_imag(doublecomplex *); + void d_cnjg(doublecomplex *, doublecomplex *); - ii = izamax_(n, &vl[ki * vl_dim1 + 1], &c__1); - i__2 = ii + ki * vl_dim1; - remax = 1. / ((d__1 = vl[i__2].r, abs(d__1)) + (d__2 = d_imag( - &vl[ii + ki * vl_dim1]), abs(d__2))); - zdscal_(n, &remax, &vl[ki * vl_dim1 + 1], &c__1); - } + /* Local variables */ + static integer i__, j, k, ii, ki, is; + static doublereal ulp; + static logical allv; + static doublereal unfl, ovfl, smin; + static logical over; + static doublereal scale; + extern logical lsame_(char *, char *); + static doublereal remax; + static logical leftv, bothv; + extern /* Subroutine */ int zgemv_(char *, integer *, integer *, + doublecomplex *, doublecomplex *, integer *, doublecomplex *, + integer *, doublecomplex *, doublecomplex *, integer *); + static logical somev; + extern /* Subroutine */ int zcopy_(integer *, doublecomplex *, integer *, + doublecomplex *, integer *), dlabad_(doublereal *, doublereal *); -/* Set back the original diagonal elements of T. */ + extern /* Subroutine */ int xerbla_(char *, integer *), zdscal_( + integer *, doublereal *, doublecomplex *, integer *); + extern integer izamax_(integer *, doublecomplex *, integer *); + static logical rightv; + extern doublereal dzasum_(integer *, doublecomplex *, integer *); + static doublereal smlnum; + extern /* Subroutine */ int zlatrs_(char *, char *, char *, char *, + integer *, doublecomplex *, integer *, doublecomplex *, + doublereal *, doublereal *, integer *); - i__2 = *n; - for (k = ki + 1; k <= i__2; ++k) { - i__3 = k + k * t_dim1; - i__4 = k + *n; - t[i__3].r = work[i__4].r, t[i__3].i = work[i__4].i; -/* L120: */ - } - ++is; -L130: - ; - } - } +/* + -- LAPACK routine (version 3.0) -- + Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., + Courant Institute, Argonne National Lab, and Rice University + June 30, 1999 - return 0; -/* End of ZTREVC */ + Purpose + ======= -} /* ztrevc_ */ + ZTREVC computes some or all of the right and/or left eigenvectors of + a complex upper triangular matrix T. + + The right eigenvector x and the left eigenvector y of T corresponding + to an eigenvalue w are defined by: + + T*x = w*x, y'*T = w*y' + + where y' denotes the conjugate transpose of the vector y. + + If all eigenvectors are requested, the routine may either return the + matrices X and/or Y of right or left eigenvectors of T, or the + products Q*X and/or Q*Y, where Q is an input unitary + matrix. If T was obtained from the Schur factorization of an + original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of + right or left eigenvectors of A. + + Arguments + ========= + + SIDE (input) CHARACTER*1 + = 'R': compute right eigenvectors only; + = 'L': compute left eigenvectors only; + = 'B': compute both right and left eigenvectors. -/* Subroutine */ int ztrti2_(char *uplo, char *diag, integer *n, - doublecomplex *a, integer *lda, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, i__1, i__2; - doublecomplex z__1; + HOWMNY (input) CHARACTER*1 + = 'A': compute all right and/or left eigenvectors; + = 'B': compute all right and/or left eigenvectors, + and backtransform them using the input matrices + supplied in VR and/or VL; + = 'S': compute selected right and/or left eigenvectors, + specified by the logical array SELECT. - /* Builtin functions */ - void z_div(doublecomplex *, doublecomplex *, doublecomplex *); + SELECT (input) LOGICAL array, dimension (N) + If HOWMNY = 'S', SELECT specifies the eigenvectors to be + computed. + If HOWMNY = 'A' or 'B', SELECT is not referenced. + To select the eigenvector corresponding to the j-th + eigenvalue, SELECT(j) must be set to .TRUE.. - /* Local variables */ - static integer j; - static doublecomplex ajj; - extern logical lsame_(char *, char *); - extern /* Subroutine */ int zscal_(integer *, doublecomplex *, - doublecomplex *, integer *); - static logical upper; - extern /* Subroutine */ int ztrmv_(char *, char *, char *, integer *, - doublecomplex *, integer *, doublecomplex *, integer *), xerbla_(char *, integer *); - static logical nounit; + N (input) INTEGER + The order of the matrix T. N >= 0. + T (input/output) COMPLEX*16 array, dimension (LDT,N) + The upper triangular matrix T. T is modified, but restored + on exit. -/* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 + LDT (input) INTEGER + The leading dimension of the array T. LDT >= max(1,N). + VL (input/output) COMPLEX*16 array, dimension (LDVL,MM) + On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must + contain an N-by-N matrix Q (usually the unitary matrix Q of + Schur vectors returned by ZHSEQR). + On exit, if SIDE = 'L' or 'B', VL contains: + if HOWMNY = 'A', the matrix Y of left eigenvectors of T; + VL is lower triangular. The i-th column + VL(i) of VL is the eigenvector corresponding + to T(i,i). + if HOWMNY = 'B', the matrix Q*Y; + if HOWMNY = 'S', the left eigenvectors of T specified by + SELECT, stored consecutively in the columns + of VL, in the same order as their + eigenvalues. + If SIDE = 'R', VL is not referenced. - Purpose - ======= + LDVL (input) INTEGER + The leading dimension of the array VL. LDVL >= max(1,N) if + SIDE = 'L' or 'B'; LDVL >= 1 otherwise. - ZTRTI2 computes the inverse of a complex upper or lower triangular - matrix. + VR (input/output) COMPLEX*16 array, dimension (LDVR,MM) + On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must + contain an N-by-N matrix Q (usually the unitary matrix Q of + Schur vectors returned by ZHSEQR). + On exit, if SIDE = 'R' or 'B', VR contains: + if HOWMNY = 'A', the matrix X of right eigenvectors of T; + VR is upper triangular. The i-th column + VR(i) of VR is the eigenvector corresponding + to T(i,i). + if HOWMNY = 'B', the matrix Q*X; + if HOWMNY = 'S', the right eigenvectors of T specified by + SELECT, stored consecutively in the columns + of VR, in the same order as their + eigenvalues. + If SIDE = 'L', VR is not referenced. - This is the Level 2 BLAS version of the algorithm. + LDVR (input) INTEGER + The leading dimension of the array VR. LDVR >= max(1,N) if + SIDE = 'R' or 'B'; LDVR >= 1 otherwise. - Arguments - ========= + MM (input) INTEGER + The number of columns in the arrays VL and/or VR. MM >= M. - UPLO (input) CHARACTER*1 - Specifies whether the matrix A is upper or lower triangular. - = 'U': Upper triangular - = 'L': Lower triangular + M (output) INTEGER + The number of columns in the arrays VL and/or VR actually + used to store the eigenvectors. If HOWMNY = 'A' or 'B', M + is set to N. Each selected eigenvector occupies one + column. - DIAG (input) CHARACTER*1 - Specifies whether or not the matrix A is unit triangular. - = 'N': Non-unit triangular - = 'U': Unit triangular + WORK (workspace) COMPLEX*16 array, dimension (2*N) - N (input) INTEGER - The order of the matrix A. N >= 0. + RWORK (workspace) DOUBLE PRECISION array, dimension (N) - A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the triangular matrix A. If UPLO = 'U', the - leading n by n upper triangular part of the array A contains - the upper triangular matrix, and the strictly lower - triangular part of A is not referenced. If UPLO = 'L', the - leading n by n lower triangular part of the array A contains - the lower triangular matrix, and the strictly upper - triangular part of A is not referenced. If DIAG = 'U', the - diagonal elements of A are also not referenced and are - assumed to be 1. + INFO (output) INTEGER + = 0: successful exit + < 0: if INFO = -i, the i-th argument had an illegal value - On exit, the (triangular) inverse of the original matrix, in - the same storage format. + Further Details + =============== - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). + The algorithm used in this program is basically backward (forward) + substitution, with scaling to make the the code robust against + possible overflow. - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -k, the k-th argument had an illegal value + Each eigenvector is normalized so that the element of largest + magnitude has magnitude 1; here the magnitude of a complex number + (x,y) is taken to be |x| + |y|. ===================================================================== - Test the input parameters. + Decode and test the input parameters */ /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; + --select; + t_dim1 = *ldt; + t_offset = 1 + t_dim1; + t -= t_offset; + vl_dim1 = *ldvl; + vl_offset = 1 + vl_dim1; + vl -= vl_offset; + vr_dim1 = *ldvr; + vr_offset = 1 + vr_dim1; + vr -= vr_offset; + --work; + --rwork; /* Function Body */ - *info = 0; - upper = lsame_(uplo, "U"); - nounit = lsame_(diag, "N"); - if (! upper && ! lsame_(uplo, "L")) { - *info = -1; - } else if (! nounit && ! lsame_(diag, "U")) { - *info = -2; - } else if (*n < 0) { - *info = -3; - } else if (*lda < max(1,*n)) { - *info = -5; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZTRTI2", &i__1); - return 0; - } + bothv = lsame_(side, "B"); + rightv = lsame_(side, "R") || bothv; + leftv = lsame_(side, "L") || bothv; - if (upper) { + allv = lsame_(howmny, "A"); + over = lsame_(howmny, "B"); + somev = lsame_(howmny, "S"); -/* Compute inverse of upper triangular matrix. */ +/* + Set M to the number of columns required to store the selected + eigenvectors. +*/ + if (somev) { + *m = 0; i__1 = *n; for (j = 1; j <= i__1; ++j) { - if (nounit) { - i__2 = j + j * a_dim1; - z_div(&z__1, &c_b60, &a[j + j * a_dim1]); - a[i__2].r = z__1.r, a[i__2].i = z__1.i; - i__2 = j + j * a_dim1; - z__1.r = -a[i__2].r, z__1.i = -a[i__2].i; - ajj.r = z__1.r, ajj.i = z__1.i; - } else { - z__1.r = -1., z__1.i = -0.; - ajj.r = z__1.r, ajj.i = z__1.i; + if (select[j]) { + ++(*m); } - -/* Compute elements 1:j-1 of j-th column. */ - - i__2 = j - 1; - ztrmv_("Upper", "No transpose", diag, &i__2, &a[a_offset], lda, & - a[j * a_dim1 + 1], &c__1); - i__2 = j - 1; - zscal_(&i__2, &ajj, &a[j * a_dim1 + 1], &c__1); /* L10: */ } } else { - -/* Compute inverse of lower triangular matrix. */ - - for (j = *n; j >= 1; --j) { - if (nounit) { - i__1 = j + j * a_dim1; - z_div(&z__1, &c_b60, &a[j + j * a_dim1]); - a[i__1].r = z__1.r, a[i__1].i = z__1.i; - i__1 = j + j * a_dim1; - z__1.r = -a[i__1].r, z__1.i = -a[i__1].i; - ajj.r = z__1.r, ajj.i = z__1.i; - } else { - z__1.r = -1., z__1.i = -0.; - ajj.r = z__1.r, ajj.i = z__1.i; - } - if (j < *n) { - -/* Compute elements j+1:n of j-th column. */ - - i__1 = *n - j; - ztrmv_("Lower", "No transpose", diag, &i__1, &a[j + 1 + (j + - 1) * a_dim1], lda, &a[j + 1 + j * a_dim1], &c__1); - i__1 = *n - j; - zscal_(&i__1, &ajj, &a[j + 1 + j * a_dim1], &c__1); - } -/* L20: */ - } + *m = *n; } - return 0; + *info = 0; + if (! rightv && ! leftv) { + *info = -1; + } else if (! allv && ! over && ! somev) { + *info = -2; + } else if (*n < 0) { + *info = -4; + } else if (*ldt < max(1,*n)) { + *info = -6; + } else if (*ldvl < 1 || leftv && *ldvl < *n) { + *info = -8; + } else if (*ldvr < 1 || rightv && *ldvr < *n) { + *info = -10; + } else if (*mm < *m) { + *info = -11; + } + if (*info != 0) { + i__1 = -(*info); + xerbla_("ZTREVC", &i__1); + return 0; + } -/* End of ZTRTI2 */ +/* Quick return if possible. */ -} /* ztrti2_ */ + if (*n == 0) { + return 0; + } -/* Subroutine */ int ztrtri_(char *uplo, char *diag, integer *n, - doublecomplex *a, integer *lda, integer *info) -{ - /* System generated locals */ - address a__1[2]; - integer a_dim1, a_offset, i__1, i__2, i__3[2], i__4, i__5; - doublecomplex z__1; - char ch__1[2]; +/* Set the constants to control overflow. */ - /* Builtin functions */ - /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen); + unfl = SAFEMINIMUM; + ovfl = 1. / unfl; + dlabad_(&unfl, &ovfl); + ulp = PRECISION; + smlnum = unfl * (*n / ulp); - /* Local variables */ - static integer j, jb, nb, nn; - extern logical lsame_(char *, char *); - static logical upper; - extern /* Subroutine */ int ztrmm_(char *, char *, char *, char *, - integer *, integer *, doublecomplex *, doublecomplex *, integer *, - doublecomplex *, integer *), - ztrsm_(char *, char *, char *, char *, integer *, integer *, - doublecomplex *, doublecomplex *, integer *, doublecomplex *, - integer *), ztrti2_(char *, char * - , integer *, doublecomplex *, integer *, integer *), xerbla_(char *, integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - static logical nounit; +/* Store the diagonal elements of T in working array WORK. */ + i__1 = *n; + for (i__ = 1; i__ <= i__1; ++i__) { + i__2 = i__ + *n; + i__3 = i__ + i__ * t_dim1; + work[i__2].r = t[i__3].r, work[i__2].i = t[i__3].i; +/* L20: */ + } /* - -- LAPACK routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 - - - Purpose - ======= + Compute 1-norm of each column of strictly upper triangular + part of T to control overflow in triangular solver. +*/ - ZTRTRI computes the inverse of a complex upper or lower triangular - matrix A. + rwork[1] = 0.; + i__1 = *n; + for (j = 2; j <= i__1; ++j) { + i__2 = j - 1; + rwork[j] = dzasum_(&i__2, &t[j * t_dim1 + 1], &c__1); +/* L30: */ + } - This is the Level 3 BLAS version of the algorithm. + if (rightv) { - Arguments - ========= +/* Compute right eigenvectors. */ - UPLO (input) CHARACTER*1 - = 'U': A is upper triangular; - = 'L': A is lower triangular. + is = *m; + for (ki = *n; ki >= 1; --ki) { - DIAG (input) CHARACTER*1 - = 'N': A is non-unit triangular; - = 'U': A is unit triangular. + if (somev) { + if (! select[ki]) { + goto L80; + } + } +/* Computing MAX */ + i__1 = ki + ki * t_dim1; + d__3 = ulp * ((d__1 = t[i__1].r, abs(d__1)) + (d__2 = d_imag(&t[ + ki + ki * t_dim1]), abs(d__2))); + smin = max(d__3,smlnum); - N (input) INTEGER - The order of the matrix A. N >= 0. + work[1].r = 1., work[1].i = 0.; - A (input/output) COMPLEX*16 array, dimension (LDA,N) - On entry, the triangular matrix A. If UPLO = 'U', the - leading N-by-N upper triangular part of the array A contains - the upper triangular matrix, and the strictly lower - triangular part of A is not referenced. If UPLO = 'L', the - leading N-by-N lower triangular part of the array A contains - the lower triangular matrix, and the strictly upper - triangular part of A is not referenced. If DIAG = 'U', the - diagonal elements of A are also not referenced and are - assumed to be 1. - On exit, the (triangular) inverse of the original matrix, in - the same storage format. +/* Form right-hand side. */ - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,N). + i__1 = ki - 1; + for (k = 1; k <= i__1; ++k) { + i__2 = k; + i__3 = k + ki * t_dim1; + z__1.r = -t[i__3].r, z__1.i = -t[i__3].i; + work[i__2].r = z__1.r, work[i__2].i = z__1.i; +/* L40: */ + } - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - > 0: if INFO = i, A(i,i) is exactly zero. The triangular - matrix is singular and its inverse can not be computed. +/* + Solve the triangular system: + (T(1:KI-1,1:KI-1) - T(KI,KI))*X = SCALE*WORK. +*/ - ===================================================================== + i__1 = ki - 1; + for (k = 1; k <= i__1; ++k) { + i__2 = k + k * t_dim1; + i__3 = k + k * t_dim1; + i__4 = ki + ki * t_dim1; + z__1.r = t[i__3].r - t[i__4].r, z__1.i = t[i__3].i - t[i__4] + .i; + t[i__2].r = z__1.r, t[i__2].i = z__1.i; + i__2 = k + k * t_dim1; + if ((d__1 = t[i__2].r, abs(d__1)) + (d__2 = d_imag(&t[k + k * + t_dim1]), abs(d__2)) < smin) { + i__3 = k + k * t_dim1; + t[i__3].r = smin, t[i__3].i = 0.; + } +/* L50: */ + } + if (ki > 1) { + i__1 = ki - 1; + zlatrs_("Upper", "No transpose", "Non-unit", "Y", &i__1, &t[ + t_offset], ldt, &work[1], &scale, &rwork[1], info); + i__1 = ki; + work[i__1].r = scale, work[i__1].i = 0.; + } - Test the input parameters. -*/ +/* Copy the vector x or Q*x to VR and normalize. */ - /* Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1; - a -= a_offset; + if (! over) { + zcopy_(&ki, &work[1], &c__1, &vr[is * vr_dim1 + 1], &c__1); - /* Function Body */ - *info = 0; - upper = lsame_(uplo, "U"); - nounit = lsame_(diag, "N"); - if (! upper && ! lsame_(uplo, "L")) { - *info = -1; - } else if (! nounit && ! lsame_(diag, "U")) { - *info = -2; - } else if (*n < 0) { - *info = -3; - } else if (*lda < max(1,*n)) { - *info = -5; - } - if (*info != 0) { - i__1 = -(*info); - xerbla_("ZTRTRI", &i__1); - return 0; - } + ii = izamax_(&ki, &vr[is * vr_dim1 + 1], &c__1); + i__1 = ii + is * vr_dim1; + remax = 1. / ((d__1 = vr[i__1].r, abs(d__1)) + (d__2 = d_imag( + &vr[ii + is * vr_dim1]), abs(d__2))); + zdscal_(&ki, &remax, &vr[is * vr_dim1 + 1], &c__1); -/* Quick return if possible */ + i__1 = *n; + for (k = ki + 1; k <= i__1; ++k) { + i__2 = k + is * vr_dim1; + vr[i__2].r = 0., vr[i__2].i = 0.; +/* L60: */ + } + } else { + if (ki > 1) { + i__1 = ki - 1; + z__1.r = scale, z__1.i = 0.; + zgemv_("N", n, &i__1, &c_b60, &vr[vr_offset], ldvr, &work[ + 1], &c__1, &z__1, &vr[ki * vr_dim1 + 1], &c__1); + } - if (*n == 0) { - return 0; - } + ii = izamax_(n, &vr[ki * vr_dim1 + 1], &c__1); + i__1 = ii + ki * vr_dim1; + remax = 1. / ((d__1 = vr[i__1].r, abs(d__1)) + (d__2 = d_imag( + &vr[ii + ki * vr_dim1]), abs(d__2))); + zdscal_(n, &remax, &vr[ki * vr_dim1 + 1], &c__1); + } -/* Check for singularity if non-unit. */ +/* Set back the original diagonal elements of T. */ - if (nounit) { - i__1 = *n; - for (*info = 1; *info <= i__1; ++(*info)) { - i__2 = *info + *info * a_dim1; - if (a[i__2].r == 0. && a[i__2].i == 0.) { - return 0; + i__1 = ki - 1; + for (k = 1; k <= i__1; ++k) { + i__2 = k + k * t_dim1; + i__3 = k + *n; + t[i__2].r = work[i__3].r, t[i__2].i = work[i__3].i; +/* L70: */ } -/* L10: */ + + --is; +L80: + ; } - *info = 0; } -/* - Determine the block size for this environment. - - Writing concatenation -*/ - i__3[0] = 1, a__1[0] = uplo; - i__3[1] = 1, a__1[1] = diag; - s_cat(ch__1, a__1, i__3, &c__2, (ftnlen)2); - nb = ilaenv_(&c__1, "ZTRTRI", ch__1, n, &c_n1, &c_n1, &c_n1, (ftnlen)6, ( - ftnlen)2); - if ((nb <= 1) || (nb >= *n)) { - -/* Use unblocked code */ + if (leftv) { - ztrti2_(uplo, diag, n, &a[a_offset], lda, info); - } else { +/* Compute left eigenvectors. */ -/* Use blocked code */ + is = 1; + i__1 = *n; + for (ki = 1; ki <= i__1; ++ki) { - if (upper) { + if (somev) { + if (! select[ki]) { + goto L130; + } + } +/* Computing MAX */ + i__2 = ki + ki * t_dim1; + d__3 = ulp * ((d__1 = t[i__2].r, abs(d__1)) + (d__2 = d_imag(&t[ + ki + ki * t_dim1]), abs(d__2))); + smin = max(d__3,smlnum); -/* Compute inverse of upper triangular matrix */ + i__2 = *n; + work[i__2].r = 1., work[i__2].i = 0.; - i__1 = *n; - i__2 = nb; - for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) { -/* Computing MIN */ - i__4 = nb, i__5 = *n - j + 1; - jb = min(i__4,i__5); +/* Form right-hand side. */ -/* Compute rows 1:j-1 of current block column */ + i__2 = *n; + for (k = ki + 1; k <= i__2; ++k) { + i__3 = k; + d_cnjg(&z__2, &t[ki + k * t_dim1]); + z__1.r = -z__2.r, z__1.i = -z__2.i; + work[i__3].r = z__1.r, work[i__3].i = z__1.i; +/* L90: */ + } - i__4 = j - 1; - ztrmm_("Left", "Upper", "No transpose", diag, &i__4, &jb, & - c_b60, &a[a_offset], lda, &a[j * a_dim1 + 1], lda); - i__4 = j - 1; - z__1.r = -1., z__1.i = -0.; - ztrsm_("Right", "Upper", "No transpose", diag, &i__4, &jb, & - z__1, &a[j + j * a_dim1], lda, &a[j * a_dim1 + 1], - lda); +/* + Solve the triangular system: + (T(KI+1:N,KI+1:N) - T(KI,KI))'*X = SCALE*WORK. +*/ -/* Compute inverse of current diagonal block */ + i__2 = *n; + for (k = ki + 1; k <= i__2; ++k) { + i__3 = k + k * t_dim1; + i__4 = k + k * t_dim1; + i__5 = ki + ki * t_dim1; + z__1.r = t[i__4].r - t[i__5].r, z__1.i = t[i__4].i - t[i__5] + .i; + t[i__3].r = z__1.r, t[i__3].i = z__1.i; + i__3 = k + k * t_dim1; + if ((d__1 = t[i__3].r, abs(d__1)) + (d__2 = d_imag(&t[k + k * + t_dim1]), abs(d__2)) < smin) { + i__4 = k + k * t_dim1; + t[i__4].r = smin, t[i__4].i = 0.; + } +/* L100: */ + } - ztrti2_("Upper", diag, &jb, &a[j + j * a_dim1], lda, info); -/* L20: */ + if (ki < *n) { + i__2 = *n - ki; + zlatrs_("Upper", "Conjugate transpose", "Non-unit", "Y", & + i__2, &t[ki + 1 + (ki + 1) * t_dim1], ldt, &work[ki + + 1], &scale, &rwork[1], info); + i__2 = ki; + work[i__2].r = scale, work[i__2].i = 0.; } - } else { -/* Compute inverse of lower triangular matrix */ +/* Copy the vector x or Q*x to VL and normalize. */ - nn = (*n - 1) / nb * nb + 1; - i__2 = -nb; - for (j = nn; i__2 < 0 ? j >= 1 : j <= 1; j += i__2) { -/* Computing MIN */ - i__1 = nb, i__4 = *n - j + 1; - jb = min(i__1,i__4); - if (j + jb <= *n) { + if (! over) { + i__2 = *n - ki + 1; + zcopy_(&i__2, &work[ki], &c__1, &vl[ki + is * vl_dim1], &c__1) + ; -/* Compute rows j+jb:n of current block column */ + i__2 = *n - ki + 1; + ii = izamax_(&i__2, &vl[ki + is * vl_dim1], &c__1) + ki - 1; + i__2 = ii + is * vl_dim1; + remax = 1. / ((d__1 = vl[i__2].r, abs(d__1)) + (d__2 = d_imag( + &vl[ii + is * vl_dim1]), abs(d__2))); + i__2 = *n - ki + 1; + zdscal_(&i__2, &remax, &vl[ki + is * vl_dim1], &c__1); - i__1 = *n - j - jb + 1; - ztrmm_("Left", "Lower", "No transpose", diag, &i__1, &jb, - &c_b60, &a[j + jb + (j + jb) * a_dim1], lda, &a[j - + jb + j * a_dim1], lda); - i__1 = *n - j - jb + 1; - z__1.r = -1., z__1.i = -0.; - ztrsm_("Right", "Lower", "No transpose", diag, &i__1, &jb, - &z__1, &a[j + j * a_dim1], lda, &a[j + jb + j * - a_dim1], lda); + i__2 = ki - 1; + for (k = 1; k <= i__2; ++k) { + i__3 = k + is * vl_dim1; + vl[i__3].r = 0., vl[i__3].i = 0.; +/* L110: */ + } + } else { + if (ki < *n) { + i__2 = *n - ki; + z__1.r = scale, z__1.i = 0.; + zgemv_("N", n, &i__2, &c_b60, &vl[(ki + 1) * vl_dim1 + 1], + ldvl, &work[ki + 1], &c__1, &z__1, &vl[ki * + vl_dim1 + 1], &c__1); } -/* Compute inverse of current diagonal block */ + ii = izamax_(n, &vl[ki * vl_dim1 + 1], &c__1); + i__2 = ii + ki * vl_dim1; + remax = 1. / ((d__1 = vl[i__2].r, abs(d__1)) + (d__2 = d_imag( + &vl[ii + ki * vl_dim1]), abs(d__2))); + zdscal_(n, &remax, &vl[ki * vl_dim1 + 1], &c__1); + } - ztrti2_("Lower", diag, &jb, &a[j + j * a_dim1], lda, info); -/* L30: */ +/* Set back the original diagonal elements of T. */ + + i__2 = *n; + for (k = ki + 1; k <= i__2; ++k) { + i__3 = k + k * t_dim1; + i__4 = k + *n; + t[i__3].r = work[i__4].r, t[i__3].i = work[i__4].i; +/* L120: */ } + + ++is; +L130: + ; } } return 0; -/* End of ZTRTRI */ +/* End of ZTREVC */ -} /* ztrtri_ */ +} /* ztrevc_ */ /* Subroutine */ int zung2r_(integer *m, integer *n, integer *k, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * @@ -23570,9 +22570,9 @@ L130: *info = 0; if (*m < 0) { *info = -1; - } else if ((*n < 0) || (*n > *m)) { + } else if (*n < 0 || *n > *m) { *info = -2; - } else if ((*k < 0) || (*k > *n)) { + } else if (*k < 0 || *k > *n) { *info = -3; } else if (*lda < max(1,*m)) { *info = -5; @@ -23776,8 +22776,8 @@ L130: *info = -1; } else if (*m < 0) { *info = -2; - } else if (((*n < 0) || (wantq && ((*n > *m) || (*n < min(*m,*k))))) || (! - wantq && ((*m > *n) || (*m < min(*n,*k))))) { + } else if (*n < 0 || wantq && (*n > *m || *n < min(*m,*k)) || ! wantq && ( + *m > *n || *m < min(*n,*k))) { *info = -3; } else if (*k < 0) { *info = -4; @@ -23809,7 +22809,7 @@ L130: /* Quick return if possible */ - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { work[1].r = 1., work[1].i = 0.; return 0; } @@ -23865,8 +22865,8 @@ L130: i__1 = *m - 1; i__2 = *m - 1; i__3 = *m - 1; - zungqr_(&i__1, &i__2, &i__3, &a[((a_dim1) << (1)) + 2], lda, & - tau[1], &work[1], lwork, &iinfo); + zungqr_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[ + 1], &work[1], lwork, &iinfo); } } } else { @@ -23920,8 +22920,8 @@ L130: i__1 = *n - 1; i__2 = *n - 1; i__3 = *n - 1; - zunglq_(&i__1, &i__2, &i__3, &a[((a_dim1) << (1)) + 2], lda, & - tau[1], &work[1], lwork, &iinfo); + zunglq_(&i__1, &i__2, &i__3, &a[(a_dim1 << 1) + 2], lda, &tau[ + 1], &work[1], lwork, &iinfo); } } } @@ -24028,9 +23028,9 @@ L130: lquery = *lwork == -1; if (*n < 0) { *info = -1; - } else if ((*ilo < 1) || (*ilo > max(1,*n))) { + } else if (*ilo < 1 || *ilo > max(1,*n)) { *info = -2; - } else if ((*ihi < min(*ilo,*n)) || (*ihi > *n)) { + } else if (*ihi < min(*ilo,*n) || *ihi > *n) { *info = -3; } else if (*lda < max(1,*n)) { *info = -5; @@ -24216,7 +23216,7 @@ L130: *info = -1; } else if (*n < *m) { *info = -2; - } else if ((*k < 0) || (*k > *m)) { + } else if (*k < 0 || *k > *m) { *info = -3; } else if (*lda < max(1,*m)) { *info = -5; @@ -24407,7 +23407,7 @@ L130: *info = -1; } else if (*n < *m) { *info = -2; - } else if ((*k < 0) || (*k > *m)) { + } else if (*k < 0 || *k > *m) { *info = -3; } else if (*lda < max(1,*m)) { *info = -5; @@ -24670,9 +23670,9 @@ L130: lquery = *lwork == -1; if (*m < 0) { *info = -1; - } else if ((*n < 0) || (*n > *m)) { + } else if (*n < 0 || *n > *m) { *info = -2; - } else if ((*k < 0) || (*k > *n)) { + } else if (*k < 0 || *k > *n) { *info = -3; } else if (*lda < max(1,*m)) { *info = -5; @@ -24966,7 +23966,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,nq)) { *info = -7; @@ -24981,11 +23981,11 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { return 0; } - if ((left && notran) || (! left && ! notran)) { + if (left && notran || ! left && ! notran) { i1 = 1; i2 = *k; i3 = 1; @@ -25181,7 +24181,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,nq)) { *info = -7; @@ -25196,11 +24196,11 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { return 0; } - if ((left && ! notran) || (! left && notran)) { + if (left && ! notran || ! left && notran) { i1 = 1; i2 = *k; i3 = 1; @@ -25450,8 +24450,7 @@ L130: } else /* if(complicated condition) */ { /* Computing MAX */ i__1 = 1, i__2 = min(nq,*k); - if ((applyq && *lda < max(1,nq)) || (! applyq && *lda < max(i__1,i__2) - )) { + if (applyq && *lda < max(1,nq) || ! applyq && *lda < max(i__1,i__2)) { *info = -8; } else if (*ldc < max(1,*m)) { *info = -11; @@ -25516,7 +24515,7 @@ L130: /* Quick return if possible */ work[1].r = 1., work[1].i = 0.; - if ((*m == 0) || (*n == 0)) { + if (*m == 0 || *n == 0) { return 0; } @@ -25580,9 +24579,9 @@ L130: i2 = 2; } i__1 = nq - 1; - zunmlq_(side, transt, &mi, &ni, &i__1, &a[((a_dim1) << (1)) + 1], - lda, &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], - lwork, &iinfo); + zunmlq_(side, transt, &mi, &ni, &i__1, &a[(a_dim1 << 1) + 1], lda, + &tau[1], &c__[i1 + i2 * c_dim1], ldc, &work[1], lwork, & + iinfo); } } work[1].r = (doublereal) lwkopt, work[1].i = 0.; @@ -25732,7 +24731,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,*k)) { *info = -7; @@ -25747,11 +24746,11 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { return 0; } - if ((left && notran) || (! left && ! notran)) { + if (left && notran || ! left && ! notran) { i1 = 1; i2 = *k; i3 = 1; @@ -25987,7 +24986,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,*k)) { *info = -7; @@ -26026,7 +25025,7 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { work[1].r = 1., work[1].i = 0.; return 0; } @@ -26052,7 +25051,7 @@ L130: iws = nw; } - if ((nb < nbmin) || (nb >= *k)) { + if (nb < nbmin || nb >= *k) { /* Use unblocked code */ @@ -26062,7 +25061,7 @@ L130: /* Use blocked code */ - if ((left && notran) || (! left && ! notran)) { + if (left && notran || ! left && ! notran) { i1 = 1; i2 = *k; i3 = nb; @@ -26297,7 +25296,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,nq)) { *info = -7; @@ -26336,7 +25335,7 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { work[1].r = 1., work[1].i = 0.; return 0; } @@ -26362,7 +25361,7 @@ L130: iws = nw; } - if ((nb < nbmin) || (nb >= *k)) { + if (nb < nbmin || nb >= *k) { /* Use unblocked code */ @@ -26372,7 +25371,7 @@ L130: /* Use blocked code */ - if ((left && notran) || (! left && ! notran)) { + if (left && notran || ! left && ! notran) { i1 = 1; i2 = *k; i3 = nb; @@ -26597,7 +25596,7 @@ L130: *info = -3; } else if (*n < 0) { *info = -4; - } else if ((*k < 0) || (*k > nq)) { + } else if (*k < 0 || *k > nq) { *info = -5; } else if (*lda < max(1,nq)) { *info = -7; @@ -26636,7 +25635,7 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (*k == 0)) { + if (*m == 0 || *n == 0 || *k == 0) { work[1].r = 1., work[1].i = 0.; return 0; } @@ -26662,7 +25661,7 @@ L130: iws = nw; } - if ((nb < nbmin) || (nb >= *k)) { + if (nb < nbmin || nb >= *k) { /* Use unblocked code */ @@ -26672,7 +25671,7 @@ L130: /* Use blocked code */ - if ((left && ! notran) || (! left && notran)) { + if (left && ! notran || ! left && notran) { i1 = 1; i2 = *k; i3 = nb; @@ -26962,7 +25961,7 @@ L130: /* Quick return if possible */ - if (((*m == 0) || (*n == 0)) || (nq == 1)) { + if (*m == 0 || *n == 0 || nq == 1) { work[1].r = 1., work[1].i = 0.; return 0; } @@ -26980,8 +25979,8 @@ L130: /* Q was determined by a call to ZHETRD with UPLO = 'U' */ i__2 = nq - 1; - zunmql_(side, trans, &mi, &ni, &i__2, &a[((a_dim1) << (1)) + 1], lda, - &tau[1], &c__[c_offset], ldc, &work[1], lwork, &iinfo); + zunmql_(side, trans, &mi, &ni, &i__2, &a[(a_dim1 << 1) + 1], lda, & + tau[1], &c__[c_offset], ldc, &work[1], lwork, &iinfo); } else { /* Q was determined by a call to ZHETRD with UPLO = 'L' */ @@ -27003,3 +26002,4 @@ L130: /* End of ZUNMTR */ } /* zunmtr_ */ + -- cgit v1.2.1