?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bump oldest cmake, compiler and CUDA versions required
[gromacs.git] / src / gromacs / linearalgebra / gmx_lapack / dbdsqr.cpp
blob7c9fdb5aae5f3e7d8469aa578147ad52ad1ddad1
1 #include <cctype>
2 #include <cmath>
3
4 #include "../gmx_blas.h"
5 #include "../gmx_lapack.h"
6
7 #include "gromacs/utility/real.h"
8
9 void
10 F77_FUNC(dbdsqr,DBDSQR)(const char *uplo,
11 int *n,
12 int *ncvt,
13 int *nru,
14 int *ncc,
15 double *d__,
16 double *e,
17 double *vt,
18 int *ldvt,
19 double *u,
20 int *ldu,
21 double *c__,
22 int *ldc,
23 double *work,
24 int *info)
25 {
26 const char xuplo = std::toupper(*uplo);
27 int c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
28 i__2;
29 double r__1, r__2, r__3, r__4;
30 double c_b15 = -.125;
31
32 int c__1 = 1;
33 double c_b49 = 1.f;
34 double c_b72 = -1.f;
35
36 double f, g, h__;
37 int i__, j, m;
38 double r__, cs;
39 int ll;
40 double sn, mu;
41 int nm1, nm12, nm13, lll;
42 double eps, sll, tol, abse;
43 int idir;
44 double abss;
45 int oldm;
46 double cosl;
47 int isub, iter;
48 double unfl, sinl, cosr, smin, smax, sinr;
49 double oldcs;
50 int oldll;
51 double shift, sigmn, oldsn = 0.;
52 int maxit;
53 double sminl;
54 double sigmx;
55 int lower;
56 double sminoa;
57 double thresh;
58 int rotate;
59 double tolmul;
60 int itmp1,itmp2;
61
62 --d__;
63 --e;
64 vt_dim1 = *ldvt;
65 vt_offset = 1 + vt_dim1;
66 vt -= vt_offset;
67 u_dim1 = *ldu;
68 u_offset = 1 + u_dim1;
69 u -= u_offset;
70 c_dim1 = *ldc;
71 c_offset = 1 + c_dim1;
72 c__ -= c_offset;
73 --work;
74
75 *info = 0;
76
77 itmp1 = (*n > 1) ? *n : 1;
78 itmp2 = (*nru > 1) ? *nru : 1;
79
80 lower = (xuplo == 'L');
81 if ( (xuplo!='U') && !lower) {
82 *info = -1;
83 } else if (*n < 0) {
84 *info = -2;
85 } else if (*ncvt < 0) {
86 *info = -3;
87 } else if (*nru < 0) {
88 *info = -4;
89 } else if (*ncc < 0) {
90 *info = -5;
91 } else if ( ((*ncvt == 0) && (*ldvt < 1)) || ((*ncvt > 0) && (*ldvt < itmp1)) ) {
92 *info = -9;
93 } else if (*ldu < itmp2) {
94 *info = -11;
95 } else if ( ((*ncc == 0) && (*ldc < 1)) || ((*ncc > 0) && (*ldc < itmp1))) {
96 *info = -13;
97 }
98 if (*info != 0) {
99 return;
100 }
101 if (*n == 0) {
102 return;
103 }
104 if (*n == 1) {
105 goto L160;
106 }
107
108 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
109
110 if (! rotate) {
111 F77_FUNC(dlasq1,DLASQ1)(n, &d__[1], &e[1], &work[1], info);
112 return;
113 }
114
115 nm1 = *n - 1;
116 nm12 = nm1 + nm1;
117 nm13 = nm12 + nm1;
118 idir = 0;
119
120 eps = GMX_DOUBLE_EPS;
121 unfl = GMX_DOUBLE_MIN/GMX_DOUBLE_EPS;
122
123 if (lower) {
124 i__1 = *n - 1;
125 for (i__ = 1; i__ <= i__1; ++i__) {
126 F77_FUNC(dlartg,DLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
127 d__[i__] = r__;
128 e[i__] = sn * d__[i__ + 1];
129 d__[i__ + 1] = cs * d__[i__ + 1];
130 work[i__] = cs;
131 work[nm1 + i__] = sn;
132 }
133
134 if (*nru > 0) {
135 F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset],
136 ldu);
137 }
138 if (*ncc > 0) {
139 F77_FUNC(dlasr,DLASR)("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset],
140 ldc);
141 }
142 }
143
144 r__3 = 100.f, r__4 = std::pow(GMX_DOUBLE_EPS,c_b15);
145 r__1 = 10.f, r__2 = (r__3<r__4) ? r__3 : r__4;
146 tolmul = (r__1>r__2) ? r__1 : r__2;
147 tol = tolmul * eps;
148 smax = 0.f;
149 i__1 = *n;
150 for (i__ = 1; i__ <= i__1; ++i__) {
151 r__2 = smax, r__3 = (r__1 = d__[i__], std::abs(r__1));
152 smax = (r__2>r__3) ? r__2 : r__3;
153 }
154 i__1 = *n - 1;
155 for (i__ = 1; i__ <= i__1; ++i__) {
156 r__2 = smax, r__3 = (r__1 = e[i__], std::abs(r__1));
157 smax = (r__2>r__3) ? r__2 : r__3;
158 }
159 sminl = 0.f;
160 if (tol >= 0.f) {
161 sminoa = std::abs(d__[1]);
162 if (sminoa == 0.f) {
163 goto L50;
164 }
165 mu = sminoa;
166 i__1 = *n;
167 for (i__ = 2; i__ <= i__1; ++i__) {
168 mu = (r__2 = d__[i__], std::abs(r__2)) * (mu / (mu + (r__1 = e[i__ -
169 1], std::abs(r__1))));
170 sminoa = (sminoa<mu) ? sminoa : mu;
171 if (sminoa == 0.f) {
172 goto L50;
173 }
174 }
175 L50:
176 sminoa /= std::sqrt((double) (*n));
177 r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl;
178 thresh = (r__1>r__2) ? r__1 : r__2;
179 } else {
180 r__1 = std::abs(tol) * smax, r__2 = *n * 6 * *n * unfl;
181 thresh = (r__1>r__2) ? r__1 : r__2;
182 }
183 maxit = *n * 6 * *n;
184 iter = 0;
185 oldll = -1;
186 oldm = -1;
187 m = *n;
188
189 L60:
190
191 if (m <= 1) {
192 goto L160;
193 }
194 if (iter > maxit) {
195 goto L200;
196 }
197
198 if (tol < 0.f && (r__1 = d__[m], std::abs(r__1)) <= thresh) {
199 d__[m] = 0.f;
200 }
201 smax = (r__1 = d__[m], std::abs(r__1));
202 smin = smax;
203 i__1 = m - 1;
204 for (lll = 1; lll <= i__1; ++lll) {
205 ll = m - lll;
206 abss = (r__1 = d__[ll], std::abs(r__1));
207 abse = (r__1 = e[ll], std::abs(r__1));
208 if (tol < 0.f && abss <= thresh) {
209 d__[ll] = 0.f;
210 }
211 if (abse <= thresh) {
212 goto L80;
213 }
214 smin = (smin<abss) ? smin : abss;
215 r__1 = (smax>abss) ? smax : abss;
216 smax = (r__1>abse) ? r__1 : abse;
217 }
218 ll = 0;
219 goto L90;
220 L80:
221 e[ll] = 0.f;
222 if (ll == m - 1) {
223 --m;
224 goto L60;
225 }
226 L90:
227 ++ll;
228 if (ll == m - 1) {
229 F77_FUNC(dlasv2,DLASV2)(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr,
230 &sinl, &cosl);
231 d__[m - 1] = sigmx;
232 e[m - 1] = 0.f;
233 d__[m] = sigmn;
234 if (*ncvt > 0) {
235 F77_FUNC(drot,DROT)(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &
236 cosr, &sinr);
237 }
238 if (*nru > 0) {
239 F77_FUNC(drot,DROT)(nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &
240 c__1, &cosl, &sinl);
241 }
242 if (*ncc > 0) {
243 F77_FUNC(drot,DROT)(ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, &
244 cosl, &sinl);
245 }
246 m += -2;
247 goto L60;
248 }
249 if (ll > oldm || m < oldll) {
250 if ((r__1 = d__[ll], std::abs(r__1)) >= (r__2 = d__[m], std::abs(r__2))) {
251 idir = 1;
252 } else {
253 idir = 2;
254 }
255 }
256 if (idir == 1) {
257
258 if( (std::abs(e[m-1]) <= std::abs(tol) * std::abs(d__[m])) ||
259 (tol<0.0 && std::abs(e[m-1])<=thresh)) {
260 e[m - 1] = 0.f;
261 goto L60;
262 }
263 if (tol >= 0.f) {
264 mu = (r__1 = d__[ll], std::abs(r__1));
265 sminl = mu;
266 i__1 = m - 1;
267 for (lll = ll; lll <= i__1; ++lll) {
268 if ((r__1 = e[lll], std::abs(r__1)) <= tol * mu) {
269 e[lll] = 0.f;
270 goto L60;
271 }
272 mu = (r__2 = d__[lll + 1], std::abs(r__2)) * (mu / (mu + (r__1 =
273 e[lll], std::abs(r__1))));
274 sminl = (sminl<mu) ? sminl : mu;
275 }
276 }
277 } else {
278 if( (std::abs(e[ll]) <= std::abs(tol)*std::abs(d__[ll])) ||
279 (tol<0.0 && std::abs(e[ll])<=thresh)) {
280 e[ll] = 0.f;
281 goto L60;
282 }
283 if (tol >= 0.f) {
284 mu = (r__1 = d__[m], std::abs(r__1));
285 sminl = mu;
286 i__1 = ll;
287 for (lll = m - 1; lll >= i__1; --lll) {
288 if ((r__1 = e[lll], std::abs(r__1)) <= tol * mu) {
289 e[lll] = 0.f;
290 goto L60;
291 }
292 mu = (r__2 = d__[lll], std::abs(r__2)) * (mu / (mu + (r__1 = e[
293 lll], std::abs(r__1))));
294 sminl = (sminl<mu) ? sminl : mu;
295 }
296 }
297 }
298 oldll = ll;
299 oldm = m;
300
301 r__1 = eps, r__2 = tol * .01f;
302 if (tol >= 0.f && *n * tol * (sminl / smax) <= ((r__1>r__2) ? r__1 : r__2)) {
303 shift = 0.f;
304 } else {
305 if (idir == 1) {
306 sll = (r__1 = d__[ll], std::abs(r__1));
307 F77_FUNC(dlas2,DLAS2)(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
308 } else {
309 sll = (r__1 = d__[m], std::abs(r__1));
310 F77_FUNC(dlas2,DLAS2)(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
311 }
312 if (sll > 0.f) {
313 r__1 = shift / sll;
314 if (r__1 * r__1 < eps) {
315 shift = 0.f;
316 }
317 }
318 }
319 iter = iter + m - ll;
320 if (shift == 0.f) {
321 if (idir == 1) {
322 cs = 1.f;
323 oldcs = 1.f;
324 i__1 = m - 1;
325 for (i__ = ll; i__ <= i__1; ++i__) {
326 r__1 = d__[i__] * cs;
327 F77_FUNC(dlartg,DLARTG)(&r__1, &e[i__], &cs, &sn, &r__);
328 if (i__ > ll) {
329 e[i__ - 1] = oldsn * r__;
330 }
331 r__1 = oldcs * r__;
332 r__2 = d__[i__ + 1] * sn;
333 F77_FUNC(dlartg,DLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
334 work[i__ - ll + 1] = cs;
335 work[i__ - ll + 1 + nm1] = sn;
336 work[i__ - ll + 1 + nm12] = oldcs;
337 work[i__ - ll + 1 + nm13] = oldsn;
338 }
339 h__ = d__[m] * cs;
340 d__[m] = h__ * oldcs;
341 e[m - 1] = h__ * oldsn;
342 if (*ncvt > 0) {
343 i__1 = m - ll + 1;
344 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
345 ll + vt_dim1], ldvt);
346 }
347 if (*nru > 0) {
348 i__1 = m - ll + 1;
349 F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
350 + 1], &u[ll * u_dim1 + 1], ldu);
351 }
352 if (*ncc > 0) {
353 i__1 = m - ll + 1;
354 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
355 + 1], &c__[ll + c_dim1], ldc);
356 }
357 if ((r__1 = e[m - 1], std::abs(r__1)) <= thresh) {
358 e[m - 1] = 0.f;
359 }
360 } else {
361 cs = 1.f;
362 oldcs = 1.f;
363 i__1 = ll + 1;
364 for (i__ = m; i__ >= i__1; --i__) {
365 r__1 = d__[i__] * cs;
366 F77_FUNC(dlartg,DLARTG)(&r__1, &e[i__ - 1], &cs, &sn, &r__);
367 if (i__ < m) {
368 e[i__] = oldsn * r__;
369 }
370 r__1 = oldcs * r__;
371 r__2 = d__[i__ - 1] * sn;
372 F77_FUNC(dlartg,DLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
373 work[i__ - ll] = cs;
374 work[i__ - ll + nm1] = -sn;
375 work[i__ - ll + nm12] = oldcs;
376 work[i__ - ll + nm13] = -oldsn;
377 }
378 h__ = d__[ll] * cs;
379 d__[ll] = h__ * oldcs;
380 e[ll] = h__ * oldsn;
381 if (*ncvt > 0) {
382 i__1 = m - ll + 1;
383 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
384 nm13 + 1], &vt[ll + vt_dim1], ldvt);
385 }
386 if (*nru > 0) {
387 i__1 = m - ll + 1;
388 F77_FUNC(dlasr,DLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
389 u_dim1 + 1], ldu);
390 }
391 if (*ncc > 0) {
392 i__1 = m - ll + 1;
393 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
394 ll + c_dim1], ldc);
395 }
396 if ((r__1 = e[ll], std::abs(r__1)) <= thresh) {
397 e[ll] = 0.f;
398 }
399 }
400 } else {
401
402 if (idir == 1) {
403 f = ((r__1 = d__[ll], std::abs(r__1)) - shift) * ( ((d__[ll] > 0) ? c_b49 : -c_b49) + shift / d__[ll]);
404 g = e[ll];
405 i__1 = m - 1;
406 for (i__ = ll; i__ <= i__1; ++i__) {
407 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosr, &sinr, &r__);
408 if (i__ > ll) {
409 e[i__ - 1] = r__;
410 }
411 f = cosr * d__[i__] + sinr * e[i__];
412 e[i__] = cosr * e[i__] - sinr * d__[i__];
413 g = sinr * d__[i__ + 1];
414 d__[i__ + 1] = cosr * d__[i__ + 1];
415 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosl, &sinl, &r__);
416 d__[i__] = r__;
417 f = cosl * e[i__] + sinl * d__[i__ + 1];
418 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
419 if (i__ < m - 1) {
420 g = sinl * e[i__ + 1];
421 e[i__ + 1] = cosl * e[i__ + 1];
422 }
423 work[i__ - ll + 1] = cosr;
424 work[i__ - ll + 1 + nm1] = sinr;
425 work[i__ - ll + 1 + nm12] = cosl;
426 work[i__ - ll + 1 + nm13] = sinl;
427 }
428 e[m - 1] = f;
429
430 if (*ncvt > 0) {
431 i__1 = m - ll + 1;
432 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
433 ll + vt_dim1], ldvt);
434 }
435 if (*nru > 0) {
436 i__1 = m - ll + 1;
437 F77_FUNC(dlasr,DLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
438 + 1], &u[ll * u_dim1 + 1], ldu);
439 }
440 if (*ncc > 0) {
441 i__1 = m - ll + 1;
442 F77_FUNC(dlasr,DLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
443 + 1], &c__[ll + c_dim1], ldc);
444 }
445 if ((r__1 = e[m - 1], std::abs(r__1)) <= thresh) {
446 e[m - 1] = 0.f;
447 }
448 } else {
449
450 f = ((r__1 = d__[m], std::abs(r__1)) - shift) * ( ((d__[m] > 0) ? c_b49 : -c_b49) + shift / d__[m]);
451 g = e[m - 1];
452 i__1 = ll + 1;
453 for (i__ = m; i__ >= i__1; --i__) {
454 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosr, &sinr, &r__);
455 if (i__ < m) {
456 e[i__] = r__;
457 }
458 f = cosr * d__[i__] + sinr * e[i__ - 1];
459 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
460 g = sinr * d__[i__ - 1];
461 d__[i__ - 1] = cosr * d__[i__ - 1];
462 F77_FUNC(dlartg,DLARTG)(&f, &g, &cosl, &sinl, &r__);
463 d__[i__] = r__;
464 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
465 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
466 if (i__ > ll + 1) {
467 g = sinl * e[i__ - 2];
468 e[i__ - 2] = cosl * e[i__ - 2];
469 }
470 work[i__ - ll] = cosr;
471 work[i__ - ll + nm1] = -sinr;
472 work[i__ - ll + nm12] = cosl;
473 work[i__ - ll + nm13] = -sinl;
474 }
475 e[ll] = f;
476
477 if ((r__1 = e[ll], std::abs(r__1)) <= thresh) {
478 e[ll] = 0.f;
479 }
480 if (*ncvt > 0) {
481 i__1 = m - ll + 1;
482 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
483 nm13 + 1], &vt[ll + vt_dim1], ldvt);
484 }
485 if (*nru > 0) {
486 i__1 = m - ll + 1;
487 F77_FUNC(dlasr,DLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
488 u_dim1 + 1], ldu);
489 }
490 if (*ncc > 0) {
491 i__1 = m - ll + 1;
492 F77_FUNC(dlasr,DLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
493 ll + c_dim1], ldc);
494 }
495 }
496 }
497
498 goto L60;
499
500 L160:
501 i__1 = *n;
502 for (i__ = 1; i__ <= i__1; ++i__) {
503 if (d__[i__] < 0.f) {
504 d__[i__] = -d__[i__];
505
506 if (*ncvt > 0) {
507 F77_FUNC(dscal,DSCAL)(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
508 }
509 }
510 }
511
512 i__1 = *n - 1;
513 for (i__ = 1; i__ <= i__1; ++i__) {
514
515 isub = 1;
516 smin = d__[1];
517 i__2 = *n + 1 - i__;
518 for (j = 2; j <= i__2; ++j) {
519 if (d__[j] <= smin) {
520 isub = j;
521 smin = d__[j];
522 }
523 }
524 if (isub != *n + 1 - i__) {
525 d__[isub] = d__[*n + 1 - i__];
526 d__[*n + 1 - i__] = smin;
527 if (*ncvt > 0) {
528 F77_FUNC(dswap,DSWAP)(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ +
529 vt_dim1], ldvt);
530 }
531 if (*nru > 0) {
532 F77_FUNC(dswap,DSWAP)(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) *
533 u_dim1 + 1], &c__1);
534 }
535 if (*ncc > 0) {
536 F77_FUNC(dswap,DSWAP)(ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ +
537 c_dim1], ldc);
538 }
539 }
540 }
541 goto L220;
542
543 L200:
544 *info = 0;
545 i__1 = *n - 1;
546 for (i__ = 1; i__ <= i__1; ++i__) {
547 if (e[i__] != 0.f) {
548 ++(*info);
549 }
550 }
551 L220:
552 return;
553
554 }
555
556
The ultimate molecular dynamics simulation package