2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2012, 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 <types/simple.h>
39 #include "gmx_lapack.h"
40 #include "lapack_limits.h"
42 void F77_FUNC(dlar1vx,DLAR1VX)(int *n,
88 for (i__ = *b1; i__ <= i__1; ++i__) {
89 if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
97 for (i__ = *bn; i__ >= i__1; --i__) {
98 if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
110 inds = (*n << 1) + 1;
117 work[inds] = lld[*b1 - 1];
119 s = work[inds] - *sigma;
121 for (i__ = *b1; i__ <= i__1; ++i__) {
122 dplus = d__[i__] + s;
123 work[i__] = ld[i__] / dplus;
124 work[inds + i__] = s * work[i__] * l[i__];
125 s = work[inds + i__] - *sigma;
128 if (! (s > 0. || s < 1.)) {
133 if (work[inds + j] > 0. || work[inds + j] < 1.) {
137 work[inds + j] = lld[j];
138 s = work[inds + j] - *sigma;
140 for (i__ = j + 1; i__ <= i__1; ++i__) {
141 dplus = d__[i__] + s;
142 work[i__] = ld[i__] / dplus;
143 if (fabs(work[i__])<GMX_DOUBLE_MIN) {
144 work[inds + i__] = lld[i__];
146 work[inds + i__] = s * work[i__] * l[i__];
148 s = work[inds + i__] - *sigma;
152 work[indp + *bn - 1] = d__[*bn] - *sigma;
154 for (i__ = *bn - 1; i__ >= i__1; --i__) {
155 dminus = lld[i__] + work[indp + i__];
156 tmp = d__[i__] / dminus;
157 work[indumn + i__] = l[i__] * tmp;
158 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
160 tmp = work[indp + r1 - 1];
161 if (! (tmp > 0. || tmp < 1.)) {
166 if (work[indp + j] > 0. || work[indp + j] < 1.) {
170 work[indp + j] = d__[j + 1] - *sigma;
172 for (i__ = j; i__ >= i__1; --i__) {
173 dminus = lld[i__] + work[indp + i__];
174 tmp = d__[i__] / dminus;
175 work[indumn + i__] = l[i__] * tmp;
176 if (fabs(tmp)<GMX_DOUBLE_MIN) {
177 work[indp + i__ - 1] = d__[i__] - *sigma;
179 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
184 *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
185 if (fabs(*mingma)<GMX_DOUBLE_MIN) {
186 *mingma = eps * work[inds + r1 - 1];
190 for (i__ = r1; i__ <= i__1; ++i__) {
191 tmp = work[inds + i__] + work[indp + i__];
192 if (fabs(tmp)<GMX_DOUBLE_MIN) {
193 tmp = eps * work[inds + i__];
195 if (fabs(tmp) < fabs(*mingma)) {
208 to = (i__1>(*b1)) ? i__1 : (*b1);
212 for (i__ = from; i__ >= i__1; --i__) {
213 z__[i__] = -(work[i__] * z__[i__ + 1]);
214 *ztz += z__[i__] * z__[i__];
216 if (fabs(z__[to]) <= eps && fabs(z__[to + 1]) <= eps) {
221 to = (i__1>*b1) ? i__1 : *b1;
227 to = (i__1<*bn) ? i__1 : *bn;
231 for (i__ = from; i__ <= i__1; ++i__) {
232 z__[i__] = -(work[indumn + i__ - 1] * z__[i__ - 1]);
233 *ztz += z__[i__] * z__[i__];
235 if (fabs(z__[to]) <= eps && fabs(z__[to - 1]) <= eps) {
240 to = (i__1<*bn) ? i__1 : *bn;
246 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
247 if (fabs(z__[i__ + 1])<GMX_DOUBLE_MIN) {
248 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
250 z__[i__] = -(work[i__] * z__[i__ + 1]);
252 if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
256 *ztz += z__[i__] * z__[i__];
260 for (i__ = *r__; i__ <= i__1; ++i__) {
261 if (fabs(z__[i__])<GMX_DOUBLE_MIN) {
262 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
264 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
266 if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
270 *ztz += z__[i__ + 1] * z__[i__ + 1];