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 "gmx_lapack.h"
40 F77_FUNC(slasd8,SLASD8)(int *icompq,
53 int difr_dim1, difr_offset, i__1, i__2;
62 float diflj, difrj, dsigj;
68 /* avoid warnings on high gcc optimization levels */
77 difr_offset = 1 + difr_dim1;
88 d__[1] = fabs(z__[1]);
92 difr[(difr_dim1 << 1) + 1] = 1.;
98 for (i__ = 1; i__ <= i__1; ++i__) {
101 /* force store and reload from memory */
102 d__2 = (*p1) + (*p2) - dsigma[i__];
111 rho = F77_FUNC(snrm2,SNRM2)(k, &z__[1], &c__1);
112 F77_FUNC(slascl,SLASCL)("G", &c__0, &c__0, &rho, &one, k, &c__1, &z__[1], k, info);
115 F77_FUNC(slaset,SLASET)("A", k, &c__1, &one, &one, &work[iwk3], k);
118 for (j = 1; j <= i__1; ++j) {
119 F77_FUNC(slasd4,SLASD4)(k, &j, &dsigma[1], &z__[1], &work[iwk1], &rho, &d__[j], &work[
125 work[iwk3i + j] = work[iwk3i + j] * work[j] * work[iwk2i + j];
127 difr[j + difr_dim1] = -work[j + 1];
129 for (i__ = 1; i__ <= i__2; ++i__) {
130 work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
131 i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
135 for (i__ = j + 1; i__ <= i__2; ++i__) {
136 work[iwk3i + i__] = work[iwk3i + i__] * work[i__] * work[iwk2i +
137 i__] / (dsigma[i__] - dsigma[j]) / (dsigma[i__] + dsigma[
143 for (i__ = 1; i__ <= i__1; ++i__) {
144 d__2 = sqrt(fabs(work[iwk3i + i__]));
145 z__[i__] = (z__[i__] > 0) ? d__2 : -d__2;
149 for (j = 1; j <= i__1; ++j) {
154 difrj = -difr[j + difr_dim1];
155 dsigjp = -dsigma[j + 1];
157 work[j] = -z__[j] / diflj / (dsigma[j] + dj);
159 for (i__ = 1; i__ <= i__2; ++i__) {
162 /* force store and reload from memory */
163 t1 = (*p1) + (*p2) - diflj;
164 work[i__] = z__[i__] / t1 / ( dsigma[i__] + dj);
167 for (i__ = j + 1; i__ <= i__2; ++i__) {
170 /* force store and reload from memory */
171 t1 = (*p1) + (*p2) - difrj;
172 work[i__] = z__[i__] / t1 / (dsigma[i__] + dj);
174 temp = F77_FUNC(snrm2,SNRM2)(k, &work[1], &c__1);
175 work[iwk2i + j] = F77_FUNC(sdot,SDOT)(k, &work[1], &c__1, &vf[1], &c__1) / temp;
176 work[iwk3i + j] = F77_FUNC(sdot,SDOT)(k, &work[1], &c__1, &vl[1], &c__1) / temp;
178 difr[j + (difr_dim1 << 1)] = temp;
182 F77_FUNC(scopy,SCOPY)(k, &work[iwk2], &c__1, &vf[1], &c__1);
183 F77_FUNC(scopy,SCOPY)(k, &work[iwk3], &c__1, &vl[1], &c__1);