3 #include "types/simple.h"
5 #include "../gmx_blas.h"
6 #include "../gmx_lapack.h"
7 #include "lapack_limits.h"
10 F77_FUNC(ssyevr,SSYEVR)(const char *jobz, const char *range, const char *uplo, int *n,
11 float *a, int *lda, float *vl, float *vu, int *
12 il, int *iu, float *abstol, int *m, float *w,
13 float *z__, int *ldz, int *isuppz, float *work,
14 int *lwork, int *iwork, int *liwork, int *info)
16 /* System generated locals */
17 int a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
23 float eps, vll, vuu, tmp1;
28 int itmp1, inddd, indee;
35 int iscale, ieeeok, indibl, indifl;
46 float posinf,neginf,negzro,newzro;
51 /* Parameter adjustments */
53 a_offset = 1 + a_dim1;
57 z_offset = 1 + z_dim1;
63 /* Check for IEEE-compliant FP */
71 negzro = fone/(neginf+fone);
77 newzro = negzro + fzero;
80 posinf = fone /newzro;
83 neginf = neginf*posinf;
86 posinf = posinf*posinf;
90 lower = (*uplo=='L' || *uplo=='l');
91 wantz = (*jobz=='V' || *jobz=='v');
92 alleig = (*range=='A' || *range=='a');
93 valeig = (*range=='V' || *range=='v');
94 indeig = (*range=='I' || *range=='i');
97 lquery = *lwork == -1 || *liwork == -1;
113 if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
115 } else if (! (alleig || valeig || indeig)) {
117 } else if (! (lower || (*uplo=='U' || *uplo=='u'))) {
121 } else if (*lda < ((*n>1) ? *n : 1) ) {
125 if (*n > 0 && *vu <= *vl) {
129 if (*il < 1 || *il > ((*n>1) ? *n : 1)) {
131 } else if (*iu < ((*n<*il) ? *n : *il) || *iu > *n) {
137 if (*ldz < 1 || (wantz && *ldz < *n)) {
139 } else if (*lwork < lwmin && ! lquery) {
141 } else if (*liwork < liwmin && ! lquery) {
149 i__1 = (nb + 1) * *n;
150 lwkopt = (i__1>lwmin) ? i__1 : lwmin;
151 work[1] = (float) lwkopt;
167 if (alleig || indeig) {
169 w[1] = a[a_dim1 + 1];
171 if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
173 w[1] = a[a_dim1 + 1];
177 z__[z_dim1 + 1] = 1.;
181 minval = GMX_FLOAT_MIN;
182 safmin = minval*(1.0+GMX_FLOAT_EPS);
185 smlnum = safmin / eps;
186 bignum = 1. / smlnum;
189 d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
190 rmax = (d__1<d__2) ? d__1 : d__2;
196 anrm = F77_FUNC(slansy,SLANSY)("M", uplo, n, &a[a_offset], lda, &work[1]);
197 if (anrm > 0. && anrm < rmin) {
200 } else if (anrm > rmax) {
207 for (j = 1; j <= i__1; ++j) {
209 F77_FUNC(sscal,SSCAL)(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
213 for (j = 1; j <= i__1; ++j) {
214 F77_FUNC(sscal,SSCAL)(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
219 abstll = *abstol * sigma;
234 llwork = *lwork - indwk + 1;
235 F77_FUNC(ssytrd,SSYTRD)(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
236 indtau], &work[indwk], &llwork, &iinfo);
239 F77_FUNC(scopy,SCOPY)(&i__1, &work[inde], &c__1, &work[indee], &c__1);
240 F77_FUNC(scopy,SCOPY)(n, &work[indd], &c__1, &work[inddd], &c__1);
242 F77_FUNC(sstegr,SSTEGR)(jobz, range, n, &work[inddd], &work[indee], vl, vu, il, iu,
243 abstol, m, &w[1], &z__[z_offset], ldz, &isuppz[1],
244 &work[indwk], lwork, &iwork[1], liwork, info);
245 if (wantz && *info == 0) {
247 llwrkn = *lwork - indwkn + 1;
248 F77_FUNC(sormtr,SORMTR)("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau]
249 , &z__[z_offset], ldz, &work[indwkn], &llwrkn, &iinfo);
262 F77_FUNC(sscal,SSCAL)(&imax, &d__1, &w[1], &c__1);
268 for (j = 1; j <= i__1; ++j) {
272 for (jj = j + 1; jj <= i__2; ++jj) {
280 itmp1 = iwork[indibl + i__ - 1];
282 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
284 iwork[indibl + j - 1] = itmp1;
285 F77_FUNC(sswap,SSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
291 work[1] = (float) lwkopt;