3 #include "types/simple.h"
5 #include "../gmx_lapack.h"
6 #include "lapack_limits.h"
8 void F77_FUNC(dlar1vx,DLAR1VX)(int *n,
54 for (i__ = *b1; i__ <= i__1; ++i__) {
55 if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
63 for (i__ = *bn; i__ >= i__1; --i__) {
64 if (*eval >= gersch[(i__ << 1) - 1] && *eval <= gersch[i__ * 2]) {
83 work[inds] = lld[*b1 - 1];
85 s = work[inds] - *sigma;
87 for (i__ = *b1; i__ <= i__1; ++i__) {
89 work[i__] = ld[i__] / dplus;
90 work[inds + i__] = s * work[i__] * l[i__];
91 s = work[inds + i__] - *sigma;
94 if (! (s > 0. || s < 1.)) {
99 if (work[inds + j] > 0. || work[inds + j] < 1.) {
103 work[inds + j] = lld[j];
104 s = work[inds + j] - *sigma;
106 for (i__ = j + 1; i__ <= i__1; ++i__) {
107 dplus = d__[i__] + s;
108 work[i__] = ld[i__] / dplus;
109 if (fabs(work[i__])<GMX_DOUBLE_MIN) {
110 work[inds + i__] = lld[i__];
112 work[inds + i__] = s * work[i__] * l[i__];
114 s = work[inds + i__] - *sigma;
118 work[indp + *bn - 1] = d__[*bn] - *sigma;
120 for (i__ = *bn - 1; i__ >= i__1; --i__) {
121 dminus = lld[i__] + work[indp + i__];
122 tmp = d__[i__] / dminus;
123 work[indumn + i__] = l[i__] * tmp;
124 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
126 tmp = work[indp + r1 - 1];
127 if (! (tmp > 0. || tmp < 1.)) {
132 if (work[indp + j] > 0. || work[indp + j] < 1.) {
136 work[indp + j] = d__[j + 1] - *sigma;
138 for (i__ = j; i__ >= i__1; --i__) {
139 dminus = lld[i__] + work[indp + i__];
140 tmp = d__[i__] / dminus;
141 work[indumn + i__] = l[i__] * tmp;
142 if (fabs(tmp)<GMX_DOUBLE_MIN) {
143 work[indp + i__ - 1] = d__[i__] - *sigma;
145 work[indp + i__ - 1] = work[indp + i__] * tmp - *sigma;
150 *mingma = work[inds + r1 - 1] + work[indp + r1 - 1];
151 if (fabs(*mingma)<GMX_DOUBLE_MIN) {
152 *mingma = eps * work[inds + r1 - 1];
156 for (i__ = r1; i__ <= i__1; ++i__) {
157 tmp = work[inds + i__] + work[indp + i__];
158 if (fabs(tmp)<GMX_DOUBLE_MIN) {
159 tmp = eps * work[inds + i__];
161 if (fabs(tmp) < fabs(*mingma)) {
174 to = (i__1>(*b1)) ? i__1 : (*b1);
178 for (i__ = from; i__ >= i__1; --i__) {
179 z__[i__] = -(work[i__] * z__[i__ + 1]);
180 *ztz += z__[i__] * z__[i__];
182 if (fabs(z__[to]) <= eps && fabs(z__[to + 1]) <= eps) {
187 to = (i__1>*b1) ? i__1 : *b1;
193 to = (i__1<*bn) ? i__1 : *bn;
197 for (i__ = from; i__ <= i__1; ++i__) {
198 z__[i__] = -(work[indumn + i__ - 1] * z__[i__ - 1]);
199 *ztz += z__[i__] * z__[i__];
201 if (fabs(z__[to]) <= eps && fabs(z__[to - 1]) <= eps) {
206 to = (i__1<*bn) ? i__1 : *bn;
212 for (i__ = *r__ - 1; i__ >= i__1; --i__) {
213 if (fabs(z__[i__ + 1])<GMX_DOUBLE_MIN) {
214 z__[i__] = -(ld[i__ + 1] / ld[i__]) * z__[i__ + 2];
216 z__[i__] = -(work[i__] * z__[i__ + 1]);
218 if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
222 *ztz += z__[i__] * z__[i__];
226 for (i__ = *r__; i__ <= i__1; ++i__) {
227 if (fabs(z__[i__])<GMX_DOUBLE_MIN) {
228 z__[i__ + 1] = -(ld[i__ - 1] / ld[i__]) * z__[i__ - 1];
230 z__[i__ + 1] = -(work[indumn + i__] * z__[i__]);
232 if (fabs(z__[i__]) <= eps && fabs(z__[i__ + 1]) <= eps) {
236 *ztz += z__[i__ + 1] * z__[i__ + 1];