1 #include <cmath>
2 #include "../gmx_blas.h"
3 #include "../gmx_lapack.h"
4 #include "lapack_limits.h"
5
6 #include "gromacs/utility/real.h"
7
8 void
F77_FUNC(dstein,DSTEIN)9 F77_FUNC(dstein,DSTEIN)(int *n,
10 double *d__,
11 double *e,
12 int *m,
13 double *w,
14 int *iblock,
15 int *isplit,
16 double *z__,
17 int *ldz,
18 double *work,
19 int *iwork,
20 int *ifail,
21 int *info)
22 {
23 int z_dim1, z_offset, i__1, i__2, i__3;
24 double d__2, d__3, d__4, d__5;
25
26 int i__, j, b1, j1, bn;
27 double xj, scl, eps, sep, nrm, tol;
28 int its;
29 double xjm, ztr, eps1;
30 int jblk, nblk;
31 int jmax;
32
33 int iseed[4], gpind, iinfo;
34 double ortol;
35 int indrv1, indrv2, indrv3, indrv4, indrv5;
36 int nrmchk;
37 int blksiz;
38 double onenrm, dtpcrt, pertol;
39 int c__2 = 2;
40 int c__1 = 1;
41 int c_n1 = -1;
42
43 --d__;
44 --e;
45 --w;
46 --iblock;
47 --isplit;
48 z_dim1 = *ldz;
49 z_offset = 1 + z_dim1;
50 z__ -= z_offset;
51 --work;
52 --iwork;
53 --ifail;
54
55 *info = 0;
56
57 xjm = 0.0;
58 i__1 = *m;
59 for (i__ = 1; i__ <= i__1; ++i__) {
60 ifail[i__] = 0;
61 }
62
63 if (*n < 0) {
64 *info = -1;
65 } else if (*m < 0 || *m > *n) {
66 *info = -4;
67 } else if (*ldz < (*n)) {
68 *info = -9;
69 } else {
70 i__1 = *m;
71 for (j = 2; j <= i__1; ++j) {
72 if (iblock[j] < iblock[j - 1]) {
73 *info = -6;
74 break;
75 }
76 if (iblock[j] == iblock[j - 1] && w[j] < w[j - 1]) {
77 *info = -5;
78 break;
79 }
80 }
81 }
82
83 if (*info != 0) {
84 return;
85 }
86
87 if (*n == 0 || *m == 0) {
88 return;
89 } else if (*n == 1) {
90 z__[z_dim1 + 1] = 1.;
91 return;
92 }
93
94 eps = GMX_DOUBLE_EPS;
95
96 for (i__ = 1; i__ <= 4; ++i__) {
97 iseed[i__ - 1] = 1;
98 }
99
100 indrv1 = 0;
101 indrv2 = indrv1 + *n;
102 indrv3 = indrv2 + *n;
103 indrv4 = indrv3 + *n;
104 indrv5 = indrv4 + *n;
105
106 j1 = 1;
107 i__1 = iblock[*m];
108 for (nblk = 1; nblk <= i__1; ++nblk) {
109
110 if (nblk == 1) {
111 b1 = 1;
112 } else {
113 b1 = isplit[nblk - 1] + 1;
114 }
115 bn = isplit[nblk];
116 blksiz = bn - b1 + 1;
117 if (blksiz == 1) {
118 continue;
119 }
120 gpind = b1;
121
122 onenrm = std::abs(d__[b1]) + std::abs(e[b1]);
123 d__3 = onenrm;
124 d__4 = std::abs(d__[bn]) + std::abs(e[bn - 1]);
125 onenrm = (d__3>d__4) ? d__3 : d__4;
126 i__2 = bn - 1;
127 for (i__ = b1 + 1; i__ <= i__2; ++i__) {
128 d__4 = onenrm;
129 d__5 = std::abs(d__[i__]) + std::abs(e[i__ - 1]) + std::abs(e[i__]);
130 onenrm = (d__4>d__5) ? d__4 : d__5;
131 }
132 ortol = onenrm * .001;
133
134 dtpcrt = std::sqrt(.1 / blksiz);
135
136 jblk = 0;
137 i__2 = *m;
138 for (j = j1; j <= i__2; ++j) {
139 if (iblock[j] != nblk) {
140 j1 = j;
141 break;
142 }
143 ++jblk;
144 xj = w[j];
145
146 if (blksiz == 1) {
147 work[indrv1 + 1] = 1.;
148 goto L120;
149 }
150
151 if (jblk > 1) {
152 eps1 = std::abs(eps * xj);
153 pertol = eps1 * 10.;
154 sep = xj - xjm;
155 if (sep < pertol) {
156 xj = xjm + pertol;
157 }
158 }
159
160 its = 0;
161 nrmchk = 0;
162
163 F77_FUNC(dlarnv,DLARNV)(&c__2, iseed, &blksiz, &work[indrv1 + 1]);
164
165 F77_FUNC(dcopy,DCOPY)(&blksiz, &d__[b1], &c__1, &work[indrv4 + 1], &c__1);
166 i__3 = blksiz - 1;
167 F77_FUNC(dcopy,DCOPY)(&i__3, &e[b1], &c__1, &work[indrv2 + 2], &c__1);
168 i__3 = blksiz - 1;
169 F77_FUNC(dcopy,DCOPY)(&i__3, &e[b1], &c__1, &work[indrv3 + 1], &c__1);
170
171 tol = 0.;
172 F77_FUNC(dlagtf,DLAGTF)(&blksiz, &work[indrv4 + 1], &xj, &work[indrv2 + 2], &work[
173 indrv3 + 1], &tol, &work[indrv5 + 1], &iwork[1], &iinfo);
174
175 L70:
176 ++its;
177 if (its > 5) {
178 goto L100;
179 }
180
181 d__2 = eps;
182 d__3 = std::abs(work[indrv4 + blksiz]);
183 scl = blksiz * onenrm * ((d__2>d__3) ? d__2 : d__3) / F77_FUNC(dasum,DASUM)(&blksiz, &work[
184 indrv1 + 1], &c__1);
185 F77_FUNC(dscal,DSCAL)(&blksiz, &scl, &work[indrv1 + 1], &c__1);
186
187 F77_FUNC(dlagts,DLAGTS)(&c_n1, &blksiz, &work[indrv4 + 1], &work[indrv2 + 2], &
188 work[indrv3 + 1], &work[indrv5 + 1], &iwork[1], &work[
189 indrv1 + 1], &tol, &iinfo);
190
191 if (jblk == 1) {
192 goto L90;
193 }
194 if (std::abs(xj - xjm) > ortol) {
195 gpind = j;
196 }
197 if (gpind != j) {
198 i__3 = j - 1;
199 for (i__ = gpind; i__ <= i__3; ++i__) {
200 ztr = -F77_FUNC(ddot,DDOT)(&blksiz, &work[indrv1 + 1], &c__1, &z__[b1 +
201 i__ * z_dim1], &c__1);
202 F77_FUNC(daxpy,DAXPY)(&blksiz, &ztr, &z__[b1 + i__ * z_dim1], &c__1, &
203 work[indrv1 + 1], &c__1);
204 }
205 }
206
207 L90:
208 jmax = F77_FUNC(idamax,IDAMAX)(&blksiz, &work[indrv1 + 1], &c__1);
209 nrm = std::abs(work[indrv1 + jmax]);
210
211 if (nrm < dtpcrt) {
212 goto L70;
213 }
214 ++nrmchk;
215 if (nrmchk < 3) {
216 goto L70;
217 }
218
219 goto L110;
220
221 L100:
222 ++(*info);
223 ifail[*info] = j;
224
225 L110:
226 scl = 1. / F77_FUNC(dnrm2,DNRM2)(&blksiz, &work[indrv1 + 1], &c__1);
227 jmax = F77_FUNC(idamax,IDAMAX)(&blksiz, &work[indrv1 + 1], &c__1);
228 if (work[indrv1 + jmax] < 0.) {
229 scl = -scl;
230 }
231 F77_FUNC(dscal,DSCAL)(&blksiz, &scl, &work[indrv1 + 1], &c__1);
232 L120:
233 i__3 = *n;
234 for (i__ = 1; i__ <= i__3; ++i__) {
235 z__[i__ + j * z_dim1] = 0.;
236 }
237 i__3 = blksiz;
238 for (i__ = 1; i__ <= i__3; ++i__) {
239 z__[b1 + i__ - 1 + j * z_dim1] = work[indrv1 + i__];
240 }
241
242 xjm = xj;
243 }
244 }
245
246 return;
247
248 }
249
250
251