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(slasd2,SLASD2)(int *nl,
67 int u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1, vt_offset;
68 int vt2_dim1, vt2_offset, i__1;
78 int psm[4], nlp1, nlp2, idxi, idxj;
89 u_offset = 1 + u_dim1;
92 vt_offset = 1 + vt_dim1;
96 u2_offset = 1 + u2_dim1;
99 vt2_offset = 1 + vt2_dim1;
115 z1 = *alpha * vt[nlp1 + nlp1 * vt_dim1];
117 for (i__ = *nl; i__ >= 1; --i__) {
118 z__[i__ + 1] = *alpha * vt[i__ + nlp1 * vt_dim1];
119 d__[i__ + 1] = d__[i__];
120 idxq[i__ + 1] = idxq[i__] + 1;
124 for (i__ = nlp2; i__ <= i__1; ++i__) {
125 z__[i__] = *beta * vt[i__ + nlp2 * vt_dim1];
129 for (i__ = 2; i__ <= i__1; ++i__) {
133 for (i__ = nlp2; i__ <= i__1; ++i__) {
138 for (i__ = nlp2; i__ <= i__1; ++i__) {
143 for (i__ = 2; i__ <= i__1; ++i__) {
144 dsigma[i__] = d__[idxq[i__]];
145 u2[i__ + u2_dim1] = z__[idxq[i__]];
146 idxc[i__] = coltyp[idxq[i__]];
149 F77_FUNC(slamrg,SLAMRG)(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
152 for (i__ = 2; i__ <= i__1; ++i__) {
154 d__[i__] = dsigma[idxi];
155 z__[i__] = u2[idxi + u2_dim1];
156 coltyp[i__] = idxc[idxi];
160 d__1 = fabs(*alpha), d__2 = fabs(*beta);
161 tol = (d__1 > d__2) ? d__1 : d__2;
163 tol = eps * 8. * ((d__2 > tol) ? d__2 : tol);
168 for (j = 2; j <= i__1; ++j) {
169 if (fabs(z__[j]) <= tol) {
189 if (fabs(z__[j]) <= tol) {
196 if (fabs(d__[j] - d__[jprev]) <= tol) {
201 tau = F77_FUNC(slapy2,SLAPY2)(&c__, &s);
207 idxjp = idxq[idx[jprev] + 1];
208 idxj = idxq[idx[j] + 1];
215 F77_FUNC(srot,SROT)(&n, &u[idxjp * u_dim1 + 1], &c__1, &u[idxj * u_dim1 + 1], &
217 F77_FUNC(srot,SROT)(&m, &vt[idxjp + vt_dim1], ldvt, &vt[idxj + vt_dim1], ldvt, &
219 if (coltyp[j] != coltyp[jprev]) {
228 u2[*k + u2_dim1] = z__[jprev];
229 dsigma[*k] = d__[jprev];
238 u2[*k + u2_dim1] = z__[jprev];
239 dsigma[*k] = d__[jprev];
244 for (j = 1; j <= 4; ++j) {
248 for (j = 2; j <= i__1; ++j) {
254 psm[1] = ctot[0] + 2;
255 psm[2] = psm[1] + ctot[1];
256 psm[3] = psm[2] + ctot[2];
259 for (j = 2; j <= i__1; ++j) {
262 idxc[psm[ct - 1]] = j;
267 for (j = 2; j <= i__1; ++j) {
270 idxj = idxq[idx[idxp[idxc[j]]] + 1];
274 F77_FUNC(scopy,SCOPY)(&n, &u[idxj * u_dim1 + 1], &c__1, &u2[j * u2_dim1 + 1], &c__1);
275 F77_FUNC(scopy,SCOPY)(&m, &vt[idxj + vt_dim1], ldvt, &vt2[j + vt2_dim1], ldvt2);
280 if (fabs(dsigma[2]) <= hlftol) {
284 z__[1] = F77_FUNC(slapy2,SLAPY2)(&z1, &z__[m]);
294 if (fabs(z1) <= tol) {
302 F77_FUNC(scopy,SCOPY)(&i__1, &u2[u2_dim1 + 2], &c__1, &z__[2], &c__1);
304 F77_FUNC(slaset,SLASET)("A", &n, &c__1, &zero, &zero, &u2[u2_offset], ldu2);
305 u2[nlp1 + u2_dim1] = 1.;
308 for (i__ = 1; i__ <= i__1; ++i__) {
309 vt[m + i__ * vt_dim1] = -s * vt[nlp1 + i__ * vt_dim1];
310 vt2[i__ * vt2_dim1 + 1] = c__ * vt[nlp1 + i__ * vt_dim1];
313 for (i__ = nlp2; i__ <= i__1; ++i__) {
314 vt2[i__ * vt2_dim1 + 1] = s * vt[m + i__ * vt_dim1];
315 vt[m + i__ * vt_dim1] = c__ * vt[m + i__ * vt_dim1];
318 F77_FUNC(scopy,SCOPY)(&m, &vt[nlp1 + vt_dim1], ldvt, &vt2[vt2_dim1 + 1], ldvt2);
321 F77_FUNC(scopy,SCOPY)(&m, &vt[m + vt_dim1], ldvt, &vt2[m + vt2_dim1], ldvt2);
326 F77_FUNC(scopy,SCOPY)(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
328 F77_FUNC(slacpy,SLACPY)("A", &n, &i__1, &u2[(*k + 1) * u2_dim1 + 1], ldu2, &u[(*k + 1)
331 F77_FUNC(slacpy,SLACPY)("A", &i__1, &m, &vt2[*k + 1 + vt2_dim1], ldvt2, &vt[*k + 1 +
334 for (j = 1; j <= 4; ++j) {
335 coltyp[j] = ctot[j - 1];