1 /* ./src_f77/dsyevr.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__10 = 10;
11 static integer c__1 = 1;
12 static integer c__2 = 2;
13 static integer c__3 = 3;
14 static integer c__4 = 4;
15 static integer c_n1 = -1;
16
dsyevr_(char * jobz,char * range,char * uplo,integer * n,doublereal * a,integer * lda,doublereal * vl,doublereal * vu,integer * il,integer * iu,doublereal * abstol,integer * m,doublereal * w,doublereal * z__,integer * ldz,integer * isuppz,doublereal * work,integer * lwork,integer * iwork,integer * liwork,integer * info,ftnlen jobz_len,ftnlen range_len,ftnlen uplo_len)17 /* Subroutine */ int dsyevr_(char *jobz, char *range, char *uplo, integer *n,
18 doublereal *a, integer *lda, doublereal *vl, doublereal *vu, integer *
19 il, integer *iu, doublereal *abstol, integer *m, doublereal *w,
20 doublereal *z__, integer *ldz, integer *isuppz, doublereal *work,
21 integer *lwork, integer *iwork, integer *liwork, integer *info,
22 ftnlen jobz_len, ftnlen range_len, ftnlen uplo_len)
23 {
24 /* System generated locals */
25 integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2;
26 doublereal d__1, d__2;
27
28 /* Builtin functions */
29 double sqrt(doublereal);
30
31 /* Local variables */
32 static integer i__, j, nb, jj;
33 static doublereal eps, vll, vuu, tmp1;
34 static integer indd, inde;
35 static doublereal anrm;
36 static integer imax;
37 static doublereal rmin, rmax;
38 static integer itmp1, inddd, indee;
39 extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *,
40 integer *);
41 static doublereal sigma;
42 extern logical lsame_(char *, char *, ftnlen, ftnlen);
43 static integer iinfo;
44 static char order[1];
45 static integer indwk;
46 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
47 doublereal *, integer *), dswap_(integer *, doublereal *, integer
48 *, doublereal *, integer *);
49 static integer lwmin;
50 static logical lower, wantz;
51 extern doublereal dlamch_(char *, ftnlen);
52 static logical alleig, indeig;
53 static integer iscale, ieeeok, indibl, indifl;
54 static logical valeig;
55 static doublereal safmin;
56 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
57 integer *, integer *, ftnlen, ftnlen);
58 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
59 static doublereal abstll, bignum;
60 static integer indtau, indisp;
61 extern /* Subroutine */ int dstein_(integer *, doublereal *, doublereal *,
62 integer *, doublereal *, integer *, integer *, doublereal *,
63 integer *, doublereal *, integer *, integer *, integer *),
64 dstegr_(char *, char *, integer *, doublereal *, doublereal *,
65 doublereal *, doublereal *, integer *, integer *, doublereal *,
66 integer *, doublereal *, doublereal *, integer *, integer *,
67 doublereal *, integer *, integer *, integer *, integer *, ftnlen,
68 ftnlen);
69 static integer indiwo, indwkn;
70 extern doublereal dlansy_(char *, char *, integer *, doublereal *,
71 integer *, doublereal *, ftnlen, ftnlen);
72 extern /* Subroutine */ int dstebz_(char *, char *, integer *, doublereal
73 *, doublereal *, integer *, integer *, doublereal *, doublereal *,
74 doublereal *, integer *, integer *, doublereal *, integer *,
75 integer *, doublereal *, integer *, integer *, ftnlen, ftnlen),
76 dsterf_(integer *, doublereal *, doublereal *, integer *);
77 static integer liwmin;
78 extern /* Subroutine */ int dormtr_(char *, char *, char *, integer *,
79 integer *, doublereal *, integer *, doublereal *, doublereal *,
80 integer *, doublereal *, integer *, integer *, ftnlen, ftnlen,
81 ftnlen);
82 static integer llwrkn, llwork, nsplit;
83 static doublereal smlnum;
84 extern /* Subroutine */ int dsytrd_(char *, integer *, doublereal *,
85 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
86 integer *, integer *, ftnlen);
87 static integer lwkopt;
88 static logical lquery;
89
90
91 /* -- LAPACK driver routine (version 3.0) -- */
92 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
93 /* Courant Institute, Argonne National Lab, and Rice University */
94 /* March 20, 2000 */
95
96 /* .. Scalar Arguments .. */
97 /* .. */
98 /* .. Array Arguments .. */
99 /* .. */
100
101 /* Purpose */
102 /* ======= */
103
104 /* DSYEVR computes selected eigenvalues and, optionally, eigenvectors */
105 /* of a real symmetric matrix T. Eigenvalues and eigenvectors can be */
106 /* selected by specifying either a range of values or a range of */
107 /* indices for the desired eigenvalues. */
108
109 /* Whenever possible, DSYEVR calls DSTEGR to compute the */
110 /* eigenspectrum using Relatively Robust Representations. DSTEGR */
111 /* computes eigenvalues by the dqds algorithm, while orthogonal */
112 /* eigenvectors are computed from various "good" L D L^T representations */
113 /* (also known as Relatively Robust Representations). Gram-Schmidt */
114 /* orthogonalization is avoided as far as possible. More specifically, */
115 /* the various steps of the algorithm are as follows. For the i-th */
116 /* unreduced block of T, */
117 /* (a) Compute T - sigma_i = L_i D_i L_i^T, such that L_i D_i L_i^T */
118 /* is a relatively robust representation, */
119 /* (b) Compute the eigenvalues, lambda_j, of L_i D_i L_i^T to high */
120 /* relative accuracy by the dqds algorithm, */
121 /* (c) If there is a cluster of close eigenvalues, "choose" sigma_i */
122 /* close to the cluster, and go to step (a), */
123 /* (d) Given the approximate eigenvalue lambda_j of L_i D_i L_i^T, */
124 /* compute the corresponding eigenvector by forming a */
125 /* rank-revealing twisted factorization. */
126 /* The desired accuracy of the output can be specified by the input */
127 /* parameter ABSTOL. */
128
129 /* For more details, see "A new O(n^2) algorithm for the symmetric */
130 /* tridiagonal eigenvalue/eigenvector problem", by Inderjit Dhillon, */
131 /* Computer Science Division Technical Report No. UCB//CSD-97-971, */
132 /* UC Berkeley, May 1997. */
133
134
135 /* Note 1 : DSYEVR calls DSTEGR when the full spectrum is requested */
136 /* on machines which conform to the ieee-754 floating point standard. */
137 /* DSYEVR calls DSTEBZ and SSTEIN on non-ieee machines and */
138 /* when partial spectrum requests are made. */
139
140 /* Normal execution of DSTEGR may create NaNs and infinities and */
141 /* hence may abort due to a floating point exception in environments */
142 /* which do not handle NaNs and infinities in the ieee standard default */
143 /* manner. */
144
145 /* Arguments */
146 /* ========= */
147
148 /* JOBZ (input) CHARACTER*1 */
149 /* = 'N': Compute eigenvalues only; */
150 /* = 'V': Compute eigenvalues and eigenvectors. */
151
152 /* RANGE (input) CHARACTER*1 */
153 /* = 'A': all eigenvalues will be found. */
154 /* = 'V': all eigenvalues in the half-open interval (VL,VU] */
155 /* will be found. */
156 /* = 'I': the IL-th through IU-th eigenvalues will be found. */
157 /* ********* For RANGE = 'V' or 'I' and IU - IL < N - 1, DSTEBZ and */
158 /* ********* DSTEIN are called */
159
160 /* UPLO (input) CHARACTER*1 */
161 /* = 'U': Upper triangle of A is stored; */
162 /* = 'L': Lower triangle of A is stored. */
163
164 /* N (input) INTEGER */
165 /* The order of the matrix A. N >= 0. */
166
167 /* A (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
168 /* On entry, the symmetric matrix A. If UPLO = 'U', the */
169 /* leading N-by-N upper triangular part of A contains the */
170 /* upper triangular part of the matrix A. If UPLO = 'L', */
171 /* the leading N-by-N lower triangular part of A contains */
172 /* the lower triangular part of the matrix A. */
173 /* On exit, the lower triangle (if UPLO='L') or the upper */
174 /* triangle (if UPLO='U') of A, including the diagonal, is */
175 /* destroyed. */
176
177 /* LDA (input) INTEGER */
178 /* The leading dimension of the array A. LDA >= max(1,N). */
179
180 /* VL (input) DOUBLE PRECISION */
181 /* VU (input) DOUBLE PRECISION */
182 /* If RANGE='V', the lower and upper bounds of the interval to */
183 /* be searched for eigenvalues. VL < VU. */
184 /* Not referenced if RANGE = 'A' or 'I'. */
185
186 /* IL (input) INTEGER */
187 /* IU (input) INTEGER */
188 /* If RANGE='I', the indices (in ascending order) of the */
189 /* smallest and largest eigenvalues to be returned. */
190 /* 1 <= IL <= IU <= N, if N > 0; IL = 1 and IU = 0 if N = 0. */
191 /* Not referenced if RANGE = 'A' or 'V'. */
192
193 /* ABSTOL (input) DOUBLE PRECISION */
194 /* The absolute error tolerance for the eigenvalues. */
195 /* An approximate eigenvalue is accepted as converged */
196 /* when it is determined to lie in an interval [a,b] */
197 /* of width less than or equal to */
198
199 /* ABSTOL + EPS * max( |a|,|b| ) , */
200
201 /* where EPS is the machine precision. If ABSTOL is less than */
202 /* or equal to zero, then EPS*|T| will be used in its place, */
203 /* where |T| is the 1-norm of the tridiagonal matrix obtained */
204 /* by reducing A to tridiagonal form. */
205
206 /* See "Computing Small Singular Values of Bidiagonal Matrices */
207 /* with Guaranteed High Relative Accuracy," by Demmel and */
208 /* Kahan, LAPACK Working Note #3. */
209
210 /* If high relative accuracy is important, set ABSTOL to */
211 /* DLAMCH( 'Safe minimum' ). Doing so will guarantee that */
212 /* eigenvalues are computed to high relative accuracy when */
213 /* possible in future releases. The current code does not */
214 /* make any guarantees about high relative accuracy, but */
215 /* furutre releases will. See J. Barlow and J. Demmel, */
216 /* "Computing Accurate Eigensystems of Scaled Diagonally */
217 /* Dominant Matrices", LAPACK Working Note #7, for a discussion */
218 /* of which matrices define their eigenvalues to high relative */
219 /* accuracy. */
220
221 /* M (output) INTEGER */
222 /* The total number of eigenvalues found. 0 <= M <= N. */
223 /* If RANGE = 'A', M = N, and if RANGE = 'I', M = IU-IL+1. */
224
225 /* W (output) DOUBLE PRECISION array, dimension (N) */
226 /* The first M elements contain the selected eigenvalues in */
227 /* ascending order. */
228
229 /* Z (output) DOUBLE PRECISION array, dimension (LDZ, max(1,M)) */
230 /* If JOBZ = 'V', then if INFO = 0, the first M columns of Z */
231 /* contain the orthonormal eigenvectors of the matrix A */
232 /* corresponding to the selected eigenvalues, with the i-th */
233 /* column of Z holding the eigenvector associated with W(i). */
234 /* If JOBZ = 'N', then Z is not referenced. */
235 /* Note: the user must ensure that at least max(1,M) columns are */
236 /* supplied in the array Z; if RANGE = 'V', the exact value of M */
237 /* is not known in advance and an upper bound must be used. */
238
239 /* LDZ (input) INTEGER */
240 /* The leading dimension of the array Z. LDZ >= 1, and if */
241 /* JOBZ = 'V', LDZ >= max(1,N). */
242
243 /* ISUPPZ (output) INTEGER array, dimension ( 2*max(1,M) ) */
244 /* The support of the eigenvectors in Z, i.e., the indices */
245 /* indicating the nonzero elements in Z. The i-th eigenvector */
246 /* is nonzero only in elements ISUPPZ( 2*i-1 ) through */
247 /* ISUPPZ( 2*i ). */
248 /* ********* Implemented only for RANGE = 'A' or 'I' and IU - IL = N - 1 */
249
250 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
251 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */
252
253 /* LWORK (input) INTEGER */
254 /* The dimension of the array WORK. LWORK >= max(1,26*N). */
255 /* For optimal efficiency, LWORK >= (NB+6)*N, */
256 /* where NB is the max of the blocksize for DSYTRD and DORMTR */
257 /* returned by ILAENV. */
258
259 /* If LWORK = -1, then a workspace query is assumed; the routine */
260 /* only calculates the optimal size of the WORK array, returns */
261 /* this value as the first entry of the WORK array, and no error */
262 /* message related to LWORK is issued by XERBLA. */
263
264 /* IWORK (workspace/output) INTEGER array, dimension (LIWORK) */
265 /* On exit, if INFO = 0, IWORK(1) returns the optimal LWORK. */
266
267 /* LIWORK (input) INTEGER */
268 /* The dimension of the array IWORK. LIWORK >= max(1,10*N). */
269
270 /* If LIWORK = -1, then a workspace query is assumed; the */
271 /* routine only calculates the optimal size of the IWORK array, */
272 /* returns this value as the first entry of the IWORK array, and */
273 /* no error message related to LIWORK is issued by XERBLA. */
274
275 /* INFO (output) INTEGER */
276 /* = 0: successful exit */
277 /* < 0: if INFO = -i, the i-th argument had an illegal value */
278 /* > 0: Internal error */
279
280 /* Further Details */
281 /* =============== */
282
283 /* Based on contributions by */
284 /* Inderjit Dhillon, IBM Almaden, USA */
285 /* Osni Marques, LBNL/NERSC, USA */
286 /* Ken Stanley, Computer Science Division, University of */
287 /* California at Berkeley, USA */
288
289 /* ===================================================================== */
290
291 /* .. Parameters .. */
292 /* .. */
293 /* .. Local Scalars .. */
294 /* .. */
295 /* .. External Functions .. */
296 /* .. */
297 /* .. External Subroutines .. */
298 /* .. */
299 /* .. Intrinsic Functions .. */
300 /* .. */
301 /* .. Executable Statements .. */
302
303 /* Test the input parameters. */
304
305 /* Parameter adjustments */
306 a_dim1 = *lda;
307 a_offset = 1 + a_dim1;
308 a -= a_offset;
309 --w;
310 z_dim1 = *ldz;
311 z_offset = 1 + z_dim1;
312 z__ -= z_offset;
313 --isuppz;
314 --work;
315 --iwork;
316
317 /* Function Body */
318 ieeeok = ilaenv_(&c__10, "DSYEVR", "N", &c__1, &c__2, &c__3, &c__4, (
319 ftnlen)6, (ftnlen)1);
320
321 lower = lsame_(uplo, "L", (ftnlen)1, (ftnlen)1);
322 wantz = lsame_(jobz, "V", (ftnlen)1, (ftnlen)1);
323 alleig = lsame_(range, "A", (ftnlen)1, (ftnlen)1);
324 valeig = lsame_(range, "V", (ftnlen)1, (ftnlen)1);
325 indeig = lsame_(range, "I", (ftnlen)1, (ftnlen)1);
326
327 lquery = *lwork == -1 || *liwork == -1;
328
329 /* Computing MAX */
330 i__1 = 1, i__2 = *n * 26;
331 lwmin = max(i__1,i__2);
332 /* Computing MAX */
333 i__1 = 1, i__2 = *n * 10;
334 liwmin = max(i__1,i__2);
335
336 *info = 0;
337 if (! (wantz || lsame_(jobz, "N", (ftnlen)1, (ftnlen)1))) {
338 *info = -1;
339 } else if (! (alleig || valeig || indeig)) {
340 *info = -2;
341 } else if (! (lower || lsame_(uplo, "U", (ftnlen)1, (ftnlen)1))) {
342 *info = -3;
343 } else if (*n < 0) {
344 *info = -4;
345 } else if (*lda < max(1,*n)) {
346 *info = -6;
347 } else {
348 if (valeig) {
349 if (*n > 0 && *vu <= *vl) {
350 *info = -8;
351 }
352 } else if (indeig) {
353 if (*il < 1 || *il > max(1,*n)) {
354 *info = -9;
355 } else if (*iu < min(*n,*il) || *iu > *n) {
356 *info = -10;
357 }
358 }
359 }
360 if (*info == 0) {
361 if (*ldz < 1 || wantz && *ldz < *n) {
362 *info = -15;
363 } else if (*lwork < lwmin && ! lquery) {
364 *info = -18;
365 } else if (*liwork < liwmin && ! lquery) {
366 *info = -20;
367 }
368 }
369
370 if (*info == 0) {
371 nb = ilaenv_(&c__1, "ZHETRD", uplo, n, &c_n1, &c_n1, &c_n1, (ftnlen)6,
372 (ftnlen)1);
373 /* Computing MAX */
374 i__1 = nb, i__2 = ilaenv_(&c__1, "ZUNMTR", uplo, n, &c_n1, &c_n1, &
375 c_n1, (ftnlen)6, (ftnlen)1);
376 nb = max(i__1,i__2);
377 /* Computing MAX */
378 i__1 = (nb + 1) * *n;
379 lwkopt = max(i__1,lwmin);
380 work[1] = (doublereal) lwkopt;
381 iwork[1] = liwmin;
382 }
383
384 if (*info != 0) {
385 i__1 = -(*info);
386 xerbla_("DSYEVR", &i__1, (ftnlen)6);
387 return 0;
388 } else if (lquery) {
389 return 0;
390 }
391
392 /* Quick return if possible */
393
394 *m = 0;
395 if (*n == 0) {
396 work[1] = 1.;
397 return 0;
398 }
399
400 if (*n == 1) {
401 work[1] = 7.;
402 if (alleig || indeig) {
403 *m = 1;
404 w[1] = a[a_dim1 + 1];
405 } else {
406 if (*vl < a[a_dim1 + 1] && *vu >= a[a_dim1 + 1]) {
407 *m = 1;
408 w[1] = a[a_dim1 + 1];
409 }
410 }
411 if (wantz) {
412 z__[z_dim1 + 1] = 1.;
413 }
414 return 0;
415 }
416
417 /* Get machine constants. */
418
419 safmin = dlamch_("Safe minimum", (ftnlen)12);
420 eps = dlamch_("Precision", (ftnlen)9);
421 smlnum = safmin / eps;
422 bignum = 1. / smlnum;
423 rmin = sqrt(smlnum);
424 /* Computing MIN */
425 d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
426 rmax = min(d__1,d__2);
427
428 /* Scale matrix to allowable range, if necessary. */
429
430 iscale = 0;
431 abstll = *abstol;
432 vll = *vl;
433 vuu = *vu;
434 anrm = dlansy_("M", uplo, n, &a[a_offset], lda, &work[1], (ftnlen)1, (
435 ftnlen)1);
436 if (anrm > 0. && anrm < rmin) {
437 iscale = 1;
438 sigma = rmin / anrm;
439 } else if (anrm > rmax) {
440 iscale = 1;
441 sigma = rmax / anrm;
442 }
443 if (iscale == 1) {
444 if (lower) {
445 i__1 = *n;
446 for (j = 1; j <= i__1; ++j) {
447 i__2 = *n - j + 1;
448 dscal_(&i__2, &sigma, &a[j + j * a_dim1], &c__1);
449 /* L10: */
450 }
451 } else {
452 i__1 = *n;
453 for (j = 1; j <= i__1; ++j) {
454 dscal_(&j, &sigma, &a[j * a_dim1 + 1], &c__1);
455 /* L20: */
456 }
457 }
458 if (*abstol > 0.) {
459 abstll = *abstol * sigma;
460 }
461 if (valeig) {
462 vll = *vl * sigma;
463 vuu = *vu * sigma;
464 }
465 }
466
467 /* Call DSYTRD to reduce symmetric matrix to tridiagonal form. */
468
469 indtau = 1;
470 inde = indtau + *n;
471 indd = inde + *n;
472 indee = indd + *n;
473 inddd = indee + *n;
474 indifl = inddd + *n;
475 indwk = indifl + *n;
476 llwork = *lwork - indwk + 1;
477 dsytrd_(uplo, n, &a[a_offset], lda, &work[indd], &work[inde], &work[
478 indtau], &work[indwk], &llwork, &iinfo, (ftnlen)1);
479
480 /* If all eigenvalues are desired */
481 /* then call DSTERF or SSTEGR and DORMTR. */
482
483 if ((alleig || indeig && *il == 1 && *iu == *n) && ieeeok == 1) {
484 if (! wantz) {
485 dcopy_(n, &work[indd], &c__1, &w[1], &c__1);
486 i__1 = *n - 1;
487 dcopy_(&i__1, &work[inde], &c__1, &work[indee], &c__1);
488 dsterf_(n, &w[1], &work[indee], info);
489 } else {
490 i__1 = *n - 1;
491 dcopy_(&i__1, &work[inde], &c__1, &work[indee], &c__1);
492 dcopy_(n, &work[indd], &c__1, &work[inddd], &c__1);
493
494 dstegr_(jobz, "A", n, &work[inddd], &work[indee], vl, vu, il, iu,
495 abstol, m, &w[1], &z__[z_offset], ldz, &isuppz[1], &work[
496 indwk], lwork, &iwork[1], liwork, info, (ftnlen)1, (
497 ftnlen)1);
498
499
500
501 /* Apply orthogonal matrix used in reduction to tridiagonal */
502 /* form to eigenvectors returned by DSTEIN. */
503
504 if (wantz && *info == 0) {
505 indwkn = inde;
506 llwrkn = *lwork - indwkn + 1;
507 dormtr_("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau]
508 , &z__[z_offset], ldz, &work[indwkn], &llwrkn, &iinfo,
509 (ftnlen)1, (ftnlen)1, (ftnlen)1);
510 }
511 }
512
513
514 if (*info == 0) {
515 *m = *n;
516 goto L30;
517 }
518 *info = 0;
519 }
520
521 /* Otherwise, call DSTEBZ and, if eigenvectors are desired, SSTEIN. */
522 /* Also call DSTEBZ and SSTEIN if SSTEGR fails. */
523
524 if (wantz) {
525 *(unsigned char *)order = 'B';
526 } else {
527 *(unsigned char *)order = 'E';
528 }
529 indifl = 1;
530 indibl = indifl + *n;
531 indisp = indibl + *n;
532 indiwo = indisp + *n;
533 dstebz_(range, order, n, &vll, &vuu, il, iu, &abstll, &work[indd], &work[
534 inde], m, &nsplit, &w[1], &iwork[indibl], &iwork[indisp], &work[
535 indwk], &iwork[indiwo], info, (ftnlen)1, (ftnlen)1);
536
537 if (wantz) {
538 dstein_(n, &work[indd], &work[inde], m, &w[1], &iwork[indibl], &iwork[
539 indisp], &z__[z_offset], ldz, &work[indwk], &iwork[indiwo], &
540 iwork[indifl], info);
541
542 /* Apply orthogonal matrix used in reduction to tridiagonal */
543 /* form to eigenvectors returned by DSTEIN. */
544
545 indwkn = inde;
546 llwrkn = *lwork - indwkn + 1;
547 dormtr_("L", uplo, "N", n, m, &a[a_offset], lda, &work[indtau], &z__[
548 z_offset], ldz, &work[indwkn], &llwrkn, &iinfo, (ftnlen)1, (
549 ftnlen)1, (ftnlen)1);
550 }
551
552 /* If matrix was scaled, then rescale eigenvalues appropriately. */
553
554 L30:
555 if (iscale == 1) {
556 if (*info == 0) {
557 imax = *m;
558 } else {
559 imax = *info - 1;
560 }
561 d__1 = 1. / sigma;
562 dscal_(&imax, &d__1, &w[1], &c__1);
563 }
564
565 /* If eigenvalues are not in order, then sort them, along with */
566 /* eigenvectors. */
567
568 if (wantz) {
569 i__1 = *m - 1;
570 for (j = 1; j <= i__1; ++j) {
571 i__ = 0;
572 tmp1 = w[j];
573 i__2 = *m;
574 for (jj = j + 1; jj <= i__2; ++jj) {
575 if (w[jj] < tmp1) {
576 i__ = jj;
577 tmp1 = w[jj];
578 }
579 /* L40: */
580 }
581
582 if (i__ != 0) {
583 itmp1 = iwork[indibl + i__ - 1];
584 w[i__] = w[j];
585 iwork[indibl + i__ - 1] = iwork[indibl + j - 1];
586 w[j] = tmp1;
587 iwork[indibl + j - 1] = itmp1;
588 dswap_(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1 + 1],
589 &c__1);
590 }
591 /* L50: */
592 }
593 }
594
595 /* Set WORK(1) to optimal workspace size. */
596
597 work[1] = (doublereal) lwkopt;
598 iwork[1] = liwmin;
599
600 return 0;
601
602 /* End of DSYEVR */
603
604 } /* dsyevr_ */
605
606