2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 #include "gmx_lapack.h"
38 #include "lapack_limits.h"
40 #include <types/simple.h>
43 F77_FUNC(dstegr,DSTEGR)(const char *jobz,
64 int z_dim1, z_offset, i__1, i__2;
70 double eps, tol, tmp, rmin, rmax;
78 int valeig,alleig,indeig;
83 int iinspl, indwrk, liwmin, nsplit;
92 z_offset = 1 + z_dim1;
98 wantz = (*jobz=='V' || *jobz=='v');
99 alleig = (*range=='A' || *range=='a');
100 valeig = (*range=='V' || *range=='v');
101 indeig = (*range=='I' || *range=='i');
103 lquery = *lwork == -1 || *liwork == -1;
108 if (! (wantz || (*jobz=='N' || *jobz=='n'))) {
110 } else if (! (alleig || valeig || indeig)) {
114 } else if (valeig && *n > 0 && *vu <= *vl) {
116 } else if (indeig && (*il < 1 || *il > *n)) {
118 } else if (indeig && (*iu < *il || *iu > *n)) {
120 } else if (*ldz < 1 || (wantz && *ldz < *n)) {
122 } else if (*lwork < lwmin && ! lquery) {
124 } else if (*liwork < liwmin && ! lquery) {
128 work[1] = (double) lwmin;
145 if (alleig || indeig) {
149 if (*vl < d__[1] && *vu >= d__[1]) {
155 z__[z_dim1 + 1] = 1.;
160 minval = GMX_DOUBLE_MIN;
161 safmin = minval*(1.0+GMX_DOUBLE_EPS);
162 eps = GMX_DOUBLE_EPS;
163 smlnum = safmin / eps;
164 bignum = 1. / smlnum;
166 d__1 = sqrt(bignum), d__2 = 1. / sqrt(sqrt(safmin));
167 rmax = (d__1<d__2) ? d__1 : d__2;
169 tnrm = F77_FUNC(dlanst,DLANST)("M", n, &d__[1], &e[1]);
170 if (tnrm > 0. && tnrm < rmin) {
172 } else if (tnrm > rmax) {
175 if ( fabs(scale-1.0)>GMX_DOUBLE_EPS) {
176 F77_FUNC(dscal,DSCAL)(n, &scale, &d__[1], &c__1);
178 F77_FUNC(dscal,DSCAL)(&i__1, &scale, &e[1], &c__1);
182 indwrk = (*n << 1) + 1;
186 iindw = (*n << 1) + 1;
190 F77_FUNC(dlarrex,DLARREX)(range, n, vl, vu, il, iu, &d__[1], &e[1], &thresh, &nsplit, &
191 iwork[iinspl], m, &w[1], &iwork[iindbl], &iwork[iindw], &work[
192 indgrs], &work[indwrk], &iwork[iindwk], &iinfo);
200 d__1 = *abstol, d__2 = (double) (*n) * eps;
201 tol = (d__1>d__2) ? d__1 : d__2;
202 F77_FUNC(dlarrvx,DLARRVX)(n, &d__[1], &e[1], &iwork[iinspl], m, &w[1], &iwork[iindbl], &
203 iwork[iindw], &work[indgrs], &tol, &z__[z_offset], ldz, &
204 isuppz[1], &work[indwrk], &iwork[iindwk], &iinfo);
212 for (j = 1; j <= i__1; ++j) {
213 itmp = iwork[iindbl + j - 1];
214 w[j] += e[iwork[iinspl + itmp - 1]];
217 if (fabs(scale-1.0)>GMX_DOUBLE_EPS) {
219 F77_FUNC(dscal,DSCAL)(m, &d__1, &w[1], &c__1);
223 for (j = 1; j <= i__1; ++j) {
227 for (jj = j + 1; jj <= i__2; ++jj) {
237 F77_FUNC(dswap,DSWAP)(n, &z__[i__ * z_dim1 + 1], &c__1, &z__[j * z_dim1
239 itmp = isuppz[(i__ << 1) - 1];
240 isuppz[(i__ << 1) - 1] = isuppz[(j << 1) - 1];
241 isuppz[(j << 1) - 1] = itmp;
242 itmp = isuppz[i__ * 2];
243 isuppz[i__ * 2] = isuppz[j * 2];
244 isuppz[j * 2] = itmp;
250 work[1] = (double) lwmin;