2 #include "gromacs/utility/real.h"
4 #include "../gmx_blas.h"
5 #include "../gmx_lapack.h"
6 #include "lapack_limits.h"
9 F77_FUNC(slasd7,SLASD7)(int *icompq,
37 int givcol_dim1, givcol_offset, givnum_dim1, givnum_offset, i__1;
44 int nlp1, nlp2, idxi, idxj;
62 givcol_dim1 = *ldgcol;
63 givcol_offset = 1 + givcol_dim1;
64 givcol -= givcol_offset;
65 givnum_dim1 = *ldgnum;
66 givnum_offset = 1 + givnum_dim1;
67 givnum -= givnum_offset;
79 z1 = *alpha * vl[nlp1];
82 for (i__ = *nl; i__ >= 1; --i__) {
83 z__[i__ + 1] = *alpha * vl[i__];
85 vf[i__ + 1] = vf[i__];
86 d__[i__ + 1] = d__[i__];
87 idxq[i__ + 1] = idxq[i__] + 1;
92 for (i__ = nlp2; i__ <= i__1; ++i__) {
93 z__[i__] = *beta * vf[i__];
97 for (i__ = nlp2; i__ <= i__1; ++i__) {
102 for (i__ = 2; i__ <= i__1; ++i__) {
103 dsigma[i__] = d__[idxq[i__]];
104 zw[i__] = z__[idxq[i__]];
105 vfw[i__] = vf[idxq[i__]];
106 vlw[i__] = vl[idxq[i__]];
109 F77_FUNC(slamrg,SLAMRG)(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
112 for (i__ = 2; i__ <= i__1; ++i__) {
114 d__[i__] = dsigma[idxi];
124 tol = (d__1>d__2) ? d__1 : d__2;
126 tol = eps * 64. * ((d__2>tol) ? d__2 : tol);
131 for (j = 2; j <= i__1; ++j) {
132 if (fabs(z__[j]) <= tol) {
151 if (fabs(z__[j]) <= tol) {
157 if (fabs(d__[j] - d__[jprev]) <= tol) {
162 tau = F77_FUNC(slapy2,SLAPY2)(c__, s);
171 idxjp = idxq[idx[jprev] + 1];
172 idxj = idxq[idx[j] + 1];
179 givcol[*givptr + (givcol_dim1 << 1)] = idxjp;
180 givcol[*givptr + givcol_dim1] = idxj;
181 givnum[*givptr + (givnum_dim1 << 1)] = *c__;
182 givnum[*givptr + givnum_dim1] = *s;
184 F77_FUNC(srot,SROT)(&c__1, &vf[jprev], &c__1, &vf[j], &c__1, c__, s);
185 F77_FUNC(srot,SROT)(&c__1, &vl[jprev], &c__1, &vl[j], &c__1, c__, s);
192 dsigma[*k] = d__[jprev];
202 dsigma[*k] = d__[jprev];
208 for (j = 2; j <= i__1; ++j) {
216 for (j = 2; j <= i__1; ++j) {
218 perm[j] = idxq[idx[jp] + 1];
219 if (perm[j] <= nlp1) {
225 F77_FUNC(scopy,SCOPY)(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
229 if (fabs(dsigma[2]) <= hlftol) {
233 z__[1] = F77_FUNC(slapy2,SLAPY2)(&z1, &z__[m]);
240 *s = -z__[m] / z__[1];
242 F77_FUNC(srot,SROT)(&c__1, &vf[m], &c__1, &vf[1], &c__1, c__, s);
243 F77_FUNC(srot,SROT)(&c__1, &vl[m], &c__1, &vl[1], &c__1, c__, s);
245 if (fabs(z1) <= tol) {
253 F77_FUNC(scopy,SCOPY)(&i__1, &zw[2], &c__1, &z__[2], &c__1);
255 F77_FUNC(scopy,SCOPY)(&i__1, &vfw[2], &c__1, &vf[2], &c__1);
257 F77_FUNC(scopy,SCOPY)(&i__1, &vlw[2], &c__1, &vl[2], &c__1);