1 #include <cmath>
2 #include "../gmx_lapack.h"
3 #include "lapack_limits.h"
4
5 #include "gromacs/utility/real.h"
6
7 void
F77_FUNC(slasq6,SLASQ6)8 F77_FUNC(slasq6,SLASQ6)(int *i0,
9 int *n0,
10 float *z__,
11 int *pp,
12 float *dmin__,
13 float *dmin1,
14 float *dmin2,
15 float *dn,
16 float *dnm1,
17 float *dnm2)
18 {
19 int i__1;
20 float d__1, d__2;
21
22 /* Local variables */
23 float d__;
24 int j4, j4p2;
25 float emin, temp;
26 const float safemin = GMX_FLOAT_MIN*(1.0+GMX_FLOAT_EPS);
27
28 --z__;
29
30 if (*n0 - *i0 - 1 <= 0) {
31 return;
32 }
33
34 j4 = (*i0 << 2) + *pp - 3;
35 emin = z__[j4 + 4];
36 d__ = z__[j4];
37 *dmin__ = d__;
38
39 if (*pp == 0) {
40 i__1 = 4*(*n0 - 3);
41 for (j4 = *i0*4; j4 <= i__1; j4 += 4) {
42 z__[j4 - 2] = d__ + z__[j4 - 1];
43 if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
44 z__[j4] = 0.;
45 d__ = z__[j4 + 1];
46 *dmin__ = d__;
47 emin = 0.;
48 } else if (safemin * z__[j4 + 1] < z__[j4 - 2] && safemin * z__[j4
49 - 2] < z__[j4 + 1]) {
50 temp = z__[j4 + 1] / z__[j4 - 2];
51 z__[j4] = z__[j4 - 1] * temp;
52 d__ *= temp;
53 } else {
54 z__[j4] = z__[j4 + 1] * (z__[j4 - 1] / z__[j4 - 2]);
55 d__ = z__[j4 + 1] * (d__ / z__[j4 - 2]);
56 }
57 if(d__<*dmin__)
58 *dmin__ = d__;
59
60 d__1 = emin, d__2 = z__[j4];
61 emin = (d__1<d__2) ? d__1 : d__2;
62 }
63 } else {
64 i__1 = 4*(*n0 - 3);
65 for (j4 = *i0 << 2; j4 <= i__1; j4 += 4) {
66 z__[j4 - 3] = d__ + z__[j4];
67 if (std::abs(z__[j4 - 3])<GMX_FLOAT_MIN) {
68 z__[j4 - 1] = 0.;
69 d__ = z__[j4 + 2];
70 *dmin__ = d__;
71 emin = 0.;
72 } else if (safemin * z__[j4 + 2] < z__[j4 - 3] && safemin * z__[j4
73 - 3] < z__[j4 + 2]) {
74 temp = z__[j4 + 2] / z__[j4 - 3];
75 z__[j4 - 1] = z__[j4] * temp;
76 d__ *= temp;
77 } else {
78 z__[j4 - 1] = z__[j4 + 2] * (z__[j4] / z__[j4 - 3]);
79 d__ = z__[j4 + 2] * (d__ / z__[j4 - 3]);
80 }
81 if(d__<*dmin__)
82 *dmin__ = d__;
83 d__1 = emin, d__2 = z__[j4 - 1];
84 emin = (d__1<d__2) ? d__1 : d__2;
85 }
86 }
87
88 *dnm2 = d__;
89 *dmin2 = *dmin__;
90 j4 = 4*(*n0 - 2) - *pp;
91 j4p2 = j4 + (*pp << 1) - 1;
92 z__[j4 - 2] = *dnm2 + z__[j4p2];
93 if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
94 z__[j4] = 0.;
95 *dnm1 = z__[j4p2 + 2];
96 *dmin__ = *dnm1;
97 emin = 0.;
98 } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
99 z__[j4p2 + 2]) {
100 temp = z__[j4p2 + 2] / z__[j4 - 2];
101 z__[j4] = z__[j4p2] * temp;
102 *dnm1 = *dnm2 * temp;
103 } else {
104 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
105 *dnm1 = z__[j4p2 + 2] * (*dnm2 / z__[j4 - 2]);
106 }
107 if(*dnm1<*dmin__)
108 *dmin__ = *dnm1;
109
110 *dmin1 = *dmin__;
111 j4 += 4;
112 j4p2 = j4 + (*pp << 1) - 1;
113 z__[j4 - 2] = *dnm1 + z__[j4p2];
114 if (std::abs(z__[j4 - 2])<GMX_FLOAT_MIN) {
115 z__[j4] = 0.;
116 *dn = z__[j4p2 + 2];
117 *dmin__ = *dn;
118 emin = 0.;
119 } else if (safemin * z__[j4p2 + 2] < z__[j4 - 2] && safemin * z__[j4 - 2] <
120 z__[j4p2 + 2]) {
121 temp = z__[j4p2 + 2] / z__[j4 - 2];
122 z__[j4] = z__[j4p2] * temp;
123 *dn = *dnm1 * temp;
124 } else {
125 z__[j4] = z__[j4p2 + 2] * (z__[j4p2] / z__[j4 - 2]);
126 *dn = z__[j4p2 + 2] * (*dnm1 / z__[j4 - 2]);
127 }
128 if(*dn<*dmin__)
129 *dmin__ = *dn;
130
131 z__[j4 + 2] = *dn;
132 z__[(*n0 << 2) - *pp] = emin;
133 return;
134
135
136 }
137