4 #include "../gmx_blas.h"
5 #include "../gmx_lapack.h"
7 #include "types/simple.h"
10 F77_FUNC(sbdsqr,SBDSQR)(const char *uplo,
26 const char xuplo = toupper(*uplo);
27 int c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1,
29 float r__1, r__2, r__3, r__4;
41 int nm1, nm12, nm13, lll;
42 float eps, sll, tol, abse;
48 float unfl, sinl, cosr, smin, smax, sinr;
51 float shift, sigmn, oldsn;
66 vt_offset = 1 + vt_dim1;
69 u_offset = 1 + u_dim1;
72 c_offset = 1 + c_dim1;
78 itmp1 = (*n > 1) ? *n : 1;
79 itmp2 = (*nru > 1) ? *nru : 1;
81 lower = (xuplo == 'L');
82 if ( (xuplo!='U') && !lower) {
86 } else if (*ncvt < 0) {
88 } else if (*nru < 0) {
90 } else if (*ncc < 0) {
92 } else if ( ((*ncvt == 0) && (*ldvt < 1)) || ((*ncvt > 0) && (*ldvt < itmp1)) ) {
94 } else if (*ldu < itmp2) {
96 } else if ( ((*ncc == 0) && (*ldc < 1)) || ((*ncc > 0) && (*ldc < itmp1))) {
109 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
112 F77_FUNC(slasq1,SLASQ1)(n, &d__[1], &e[1], &work[1], info);
122 unfl = GMX_FLOAT_MIN/GMX_FLOAT_EPS;
126 for (i__ = 1; i__ <= i__1; ++i__) {
127 F77_FUNC(slartg,SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
129 e[i__] = sn * d__[i__ + 1];
130 d__[i__ + 1] = cs * d__[i__ + 1];
132 work[nm1 + i__] = sn;
136 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset],
140 F77_FUNC(slasr,SLASR)("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset],
145 r__3 = 100.f, r__4 = pow(GMX_FLOAT_EPS,c_b15);
146 r__1 = 10.f, r__2 = (r__3<r__4) ? r__3 : r__4;
147 tolmul = (r__1>r__2) ? r__1 : r__2;
151 for (i__ = 1; i__ <= i__1; ++i__) {
152 r__2 = smax, r__3 = (r__1 = d__[i__], fabs(r__1));
153 smax = (r__2>r__3) ? r__2 : r__3;
156 for (i__ = 1; i__ <= i__1; ++i__) {
157 r__2 = smax, r__3 = (r__1 = e[i__], fabs(r__1));
158 smax = (r__2>r__3) ? r__2 : r__3;
162 sminoa = fabs(d__[1]);
168 for (i__ = 2; i__ <= i__1; ++i__) {
169 mu = (r__2 = d__[i__], fabs(r__2)) * (mu / (mu + (r__1 = e[i__ -
171 sminoa = (sminoa<mu) ? sminoa : mu;
177 sminoa /= sqrt((float) (*n));
178 r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl;
179 thresh = (r__1>r__2) ? r__1 : r__2;
181 r__1 = fabs(tol) * smax, r__2 = *n * 6 * *n * unfl;
182 thresh = (r__1>r__2) ? r__1 : r__2;
199 if (tol < 0.f && (r__1 = d__[m], fabs(r__1)) <= thresh) {
202 smax = (r__1 = d__[m], fabs(r__1));
205 for (lll = 1; lll <= i__1; ++lll) {
207 abss = (r__1 = d__[ll], fabs(r__1));
208 abse = (r__1 = e[ll], fabs(r__1));
209 if (tol < 0.f && abss <= thresh) {
212 if (abse <= thresh) {
215 smin = (smin<abss) ? smin : abss;
216 r__1 = (smax>abss) ? smax : abss;
217 smax = (r__1>abse) ? r__1 : abse;
230 F77_FUNC(slasv2,SLASV2)(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr,
236 F77_FUNC(srot,SROT)(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &
240 F77_FUNC(srot,SROT)(nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &
244 F77_FUNC(srot,SROT)(ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, &
250 if (ll > oldm || m < oldll) {
251 if ((r__1 = d__[ll], fabs(r__1)) >= (r__2 = d__[m], fabs(r__2))) {
259 if( (fabs(e[m-1]) <= fabs(tol) * fabs(d__[m])) ||
260 (tol<0.0 && fabs(e[m-1])<=thresh)) {
265 mu = (r__1 = d__[ll], fabs(r__1));
268 for (lll = ll; lll <= i__1; ++lll) {
269 if ((r__1 = e[lll], fabs(r__1)) <= tol * mu) {
274 mu = (r__2 = d__[lll + 1], fabs(r__2)) * (mu / (mu + (r__1 =
275 e[lll], fabs(r__1))));
276 sminl = (sminl<mu) ? sminl : mu;
280 if( (fabs(e[ll]) <= fabs(tol)*fabs(d__[ll])) ||
281 (tol<0.0 && fabs(e[ll])<=thresh)) {
286 mu = (r__1 = d__[m], fabs(r__1));
289 for (lll = m - 1; lll >= i__1; --lll) {
290 if ((r__1 = e[lll], fabs(r__1)) <= tol * mu) {
295 mu = (r__2 = d__[lll], fabs(r__2)) * (mu / (mu + (r__1 = e[
297 sminl = (sminl<mu) ? sminl : mu;
304 r__1 = eps, r__2 = tol * .01f;
305 if (tol >= 0.f && *n * tol * (sminl / smax) <= ((r__1>r__2) ? r__1 : r__2)) {
309 sll = (r__1 = d__[ll], fabs(r__1));
310 F77_FUNC(slas2,SLAS2)(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
312 sll = (r__1 = d__[m], fabs(r__1));
313 F77_FUNC(slas2,SLAS2)(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
317 if (r__1 * r__1 < eps) {
322 iter = iter + m - ll;
328 for (i__ = ll; i__ <= i__1; ++i__) {
329 r__1 = d__[i__] * cs;
330 F77_FUNC(slartg,SLARTG)(&r__1, &e[i__], &cs, &sn, &r__);
332 e[i__ - 1] = oldsn * r__;
335 r__2 = d__[i__ + 1] * sn;
336 F77_FUNC(slartg,SLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
337 work[i__ - ll + 1] = cs;
338 work[i__ - ll + 1 + nm1] = sn;
339 work[i__ - ll + 1 + nm12] = oldcs;
340 work[i__ - ll + 1 + nm13] = oldsn;
343 d__[m] = h__ * oldcs;
344 e[m - 1] = h__ * oldsn;
347 F77_FUNC(slasr,SLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
348 ll + vt_dim1], ldvt);
352 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
353 + 1], &u[ll * u_dim1 + 1], ldu);
357 F77_FUNC(slasr,SLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
358 + 1], &c__[ll + c_dim1], ldc);
360 if ((r__1 = e[m - 1], fabs(r__1)) <= thresh) {
367 for (i__ = m; i__ >= i__1; --i__) {
368 r__1 = d__[i__] * cs;
369 F77_FUNC(slartg,SLARTG)(&r__1, &e[i__ - 1], &cs, &sn, &r__);
371 e[i__] = oldsn * r__;
374 r__2 = d__[i__ - 1] * sn;
375 F77_FUNC(slartg,SLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
377 work[i__ - ll + nm1] = -sn;
378 work[i__ - ll + nm12] = oldcs;
379 work[i__ - ll + nm13] = -oldsn;
382 d__[ll] = h__ * oldcs;
386 F77_FUNC(slasr,SLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
387 nm13 + 1], &vt[ll + vt_dim1], ldvt);
391 F77_FUNC(slasr,SLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
396 F77_FUNC(slasr,SLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
399 if ((r__1 = e[ll], fabs(r__1)) <= thresh) {
406 f = ((r__1 = d__[ll], fabs(r__1)) - shift) * ( ((d__[ll] > 0) ? c_b49 : -c_b49) + shift / d__[ll]);
409 for (i__ = ll; i__ <= i__1; ++i__) {
410 F77_FUNC(slartg,SLARTG)(&f, &g, &cosr, &sinr, &r__);
414 f = cosr * d__[i__] + sinr * e[i__];
415 e[i__] = cosr * e[i__] - sinr * d__[i__];
416 g = sinr * d__[i__ + 1];
417 d__[i__ + 1] = cosr * d__[i__ + 1];
418 F77_FUNC(slartg,SLARTG)(&f, &g, &cosl, &sinl, &r__);
420 f = cosl * e[i__] + sinl * d__[i__ + 1];
421 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
423 g = sinl * e[i__ + 1];
424 e[i__ + 1] = cosl * e[i__ + 1];
426 work[i__ - ll + 1] = cosr;
427 work[i__ - ll + 1 + nm1] = sinr;
428 work[i__ - ll + 1 + nm12] = cosl;
429 work[i__ - ll + 1 + nm13] = sinl;
435 F77_FUNC(slasr,SLASR)("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[
436 ll + vt_dim1], ldvt);
440 F77_FUNC(slasr,SLASR)("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13
441 + 1], &u[ll * u_dim1 + 1], ldu);
445 F77_FUNC(slasr,SLASR)("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13
446 + 1], &c__[ll + c_dim1], ldc);
448 if ((r__1 = e[m - 1], fabs(r__1)) <= thresh) {
453 f = ((r__1 = d__[m], fabs(r__1)) - shift) * ( ((d__[m] > 0) ? c_b49 : -c_b49) + shift / d__[m]);
456 for (i__ = m; i__ >= i__1; --i__) {
457 F77_FUNC(slartg,SLARTG)(&f, &g, &cosr, &sinr, &r__);
461 f = cosr * d__[i__] + sinr * e[i__ - 1];
462 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
463 g = sinr * d__[i__ - 1];
464 d__[i__ - 1] = cosr * d__[i__ - 1];
465 F77_FUNC(slartg,SLARTG)(&f, &g, &cosl, &sinl, &r__);
467 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
468 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
470 g = sinl * e[i__ - 2];
471 e[i__ - 2] = cosl * e[i__ - 2];
473 work[i__ - ll] = cosr;
474 work[i__ - ll + nm1] = -sinr;
475 work[i__ - ll + nm12] = cosl;
476 work[i__ - ll + nm13] = -sinl;
480 if ((r__1 = e[ll], fabs(r__1)) <= thresh) {
485 F77_FUNC(slasr,SLASR)("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[
486 nm13 + 1], &vt[ll + vt_dim1], ldvt);
490 F77_FUNC(slasr,SLASR)("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll *
495 F77_FUNC(slasr,SLASR)("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[
505 for (i__ = 1; i__ <= i__1; ++i__) {
506 if (d__[i__] < 0.f) {
507 d__[i__] = -d__[i__];
510 F77_FUNC(sscal,SSCAL)(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
516 for (i__ = 1; i__ <= i__1; ++i__) {
521 for (j = 2; j <= i__2; ++j) {
522 if (d__[j] <= smin) {
527 if (isub != *n + 1 - i__) {
528 d__[isub] = d__[*n + 1 - i__];
529 d__[*n + 1 - i__] = smin;
531 F77_FUNC(sswap,SSWAP)(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ +
535 F77_FUNC(sswap,SSWAP)(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) *
539 F77_FUNC(sswap,SSWAP)(ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ +
549 for (i__ = 1; i__ <= i__1; ++i__) {