1 /* ../netlib/dbdsqr.f -- translated by f2c (version 20100827). You must link the resulting object file with libf2c: on Microsoft Windows system, link with libf2c.lib;
2 on Linux or Unix systems, link with .../path/to/libf2c.a -lm or, if you install libf2c.a in a standard place, with -lf2c -lm -- in that order, at the end of the command line, as in cc *.o -lf2c -lm Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., http://www.netlib.org/f2c/libf2c.zip */
3 #include "FLA_f2c.h" /* Table of constant values */
4 static doublereal c_b15 = -.125;
5 static integer c__1 = 1;
6 static doublereal c_b49 = 1.;
7 static doublereal c_b72 = -1.;
8 /* > \brief \b DBDSQR */
9 /* =========== DOCUMENTATION =========== */
10 /* Online html documentation available at */
11 /* http://www.netlib.org/lapack/explore-html/ */
12 /* > \htmlonly */
13 /* > Download DBDSQR + dependencies */
14 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dbdsqr. f"> */
15 /* > [TGZ]</a> */
16 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dbdsqr. f"> */
17 /* > [ZIP]</a> */
18 /* > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dbdsqr. f"> */
19 /* > [TXT]</a> */
20 /* > \endhtmlonly */
21 /* Definition: */
22 /* =========== */
23 /* SUBROUTINE DBDSQR( UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, */
24 /* LDU, C, LDC, WORK, INFO ) */
25 /* .. Scalar Arguments .. */
26 /* CHARACTER UPLO */
27 /* INTEGER INFO, LDC, LDU, LDVT, N, NCC, NCVT, NRU */
28 /* .. */
29 /* .. Array Arguments .. */
30 /* DOUBLE PRECISION C( LDC, * ), D( * ), E( * ), U( LDU, * ), */
31 /* $ VT( LDVT, * ), WORK( * ) */
32 /* .. */
33 /* > \par Purpose: */
34 /* ============= */
35 /* > */
36 /* > \verbatim */
37 /* > */
38 /* > DBDSQR computes the singular values and, optionally, the right and/or */
39 /* > left singular vectors from the singular value decomposition (SVD) of */
40 /* > a real N-by-N (upper or lower) bidiagonal matrix B using the implicit */
41 /* > zero-shift QR algorithm. The SVD of B has the form */
42 /* > */
43 /* > B = Q * S * P**T */
44 /* > */
45 /* > where S is the diagonal matrix of singular values, Q is an orthogonal */
46 /* > matrix of left singular vectors, and P is an orthogonal matrix of */
47 /* > right singular vectors. If left singular vectors are requested, this */
48 /* > subroutine actually returns U*Q instead of Q, and, if right singular */
49 /* > vectors are requested, this subroutine returns P**T*VT instead of */
50 /* > P**T, for given real input matrices U and VT. When U and VT are the */
51 /* > orthogonal matrices that reduce a general matrix A to bidiagonal */
52 /* > form: A = U*B*VT, as computed by DGEBRD, then */
53 /* > */
54 /* > A = (U*Q) * S * (P**T*VT) */
55 /* > */
56 /* > is the SVD of A. Optionally, the subroutine may also compute Q**T*C */
57 /* > for a given real input matrix C. */
58 /* > */
59 /* > See "Computing Small Singular Values of Bidiagonal Matrices With */
60 /* > Guaranteed High Relative Accuracy," by J. Demmel and W. Kahan, */
61 /* > LAPACK Working Note #3 (or SIAM J. Sci. Statist. Comput. vol. 11, */
62 /* > no. 5, pp. 873-912, Sept 1990) and */
63 /* > "Accurate singular values and differential qd algorithms," by */
64 /* > B. Parlett and V. Fernando, Technical Report CPAM-554, Mathematics */
65 /* > Department, University of California at Berkeley, July 1992 */
66 /* > for a detailed description of the algorithm. */
67 /* > \endverbatim */
68 /* Arguments: */
69 /* ========== */
70 /* > \param[in] UPLO */
71 /* > \verbatim */
72 /* > UPLO is CHARACTER*1 */
73 /* > = 'U': B is upper bidiagonal;
74 */
75 /* > = 'L': B is lower bidiagonal. */
76 /* > \endverbatim */
77 /* > */
78 /* > \param[in] N */
79 /* > \verbatim */
80 /* > N is INTEGER */
81 /* > The order of the matrix B. N >= 0. */
82 /* > \endverbatim */
83 /* > */
84 /* > \param[in] NCVT */
85 /* > \verbatim */
86 /* > NCVT is INTEGER */
87 /* > The number of columns of the matrix VT. NCVT >= 0. */
88 /* > \endverbatim */
89 /* > */
90 /* > \param[in] NRU */
91 /* > \verbatim */
92 /* > NRU is INTEGER */
93 /* > The number of rows of the matrix U. NRU >= 0. */
94 /* > \endverbatim */
95 /* > */
96 /* > \param[in] NCC */
97 /* > \verbatim */
98 /* > NCC is INTEGER */
99 /* > The number of columns of the matrix C. NCC >= 0. */
100 /* > \endverbatim */
101 /* > */
102 /* > \param[in,out] D */
103 /* > \verbatim */
104 /* > D is DOUBLE PRECISION array, dimension (N) */
105 /* > On entry, the n diagonal elements of the bidiagonal matrix B. */
106 /* > On exit, if INFO=0, the singular values of B in decreasing */
107 /* > order. */
108 /* > \endverbatim */
109 /* > */
110 /* > \param[in,out] E */
111 /* > \verbatim */
112 /* > E is DOUBLE PRECISION array, dimension (N-1) */
113 /* > On entry, the N-1 offdiagonal elements of the bidiagonal */
114 /* > matrix B. */
115 /* > On exit, if INFO = 0, E is destroyed;
116 if INFO > 0, D and E */
117 /* > will contain the diagonal and superdiagonal elements of a */
118 /* > bidiagonal matrix orthogonally equivalent to the one given */
119 /* > as input. */
120 /* > \endverbatim */
121 /* > */
122 /* > \param[in,out] VT */
123 /* > \verbatim */
124 /* > VT is DOUBLE PRECISION array, dimension (LDVT, NCVT) */
125 /* > On entry, an N-by-NCVT matrix VT. */
126 /* > On exit, VT is overwritten by P**T * VT. */
127 /* > Not referenced if NCVT = 0. */
128 /* > \endverbatim */
129 /* > */
130 /* > \param[in] LDVT */
131 /* > \verbatim */
132 /* > LDVT is INTEGER */
133 /* > The leading dimension of the array VT. */
134 /* > LDVT >= max(1,N) if NCVT > 0;
135 LDVT >= 1 if NCVT = 0. */
136 /* > \endverbatim */
137 /* > */
138 /* > \param[in,out] U */
139 /* > \verbatim */
140 /* > U is DOUBLE PRECISION array, dimension (LDU, N) */
141 /* > On entry, an NRU-by-N matrix U. */
142 /* > On exit, U is overwritten by U * Q. */
143 /* > Not referenced if NRU = 0. */
144 /* > \endverbatim */
145 /* > */
146 /* > \param[in] LDU */
147 /* > \verbatim */
148 /* > LDU is INTEGER */
149 /* > The leading dimension of the array U. LDU >= max(1,NRU). */
150 /* > \endverbatim */
151 /* > */
152 /* > \param[in,out] C */
153 /* > \verbatim */
154 /* > C is DOUBLE PRECISION array, dimension (LDC, NCC) */
155 /* > On entry, an N-by-NCC matrix C. */
156 /* > On exit, C is overwritten by Q**T * C. */
157 /* > Not referenced if NCC = 0. */
158 /* > \endverbatim */
159 /* > */
160 /* > \param[in] LDC */
161 /* > \verbatim */
162 /* > LDC is INTEGER */
163 /* > The leading dimension of the array C. */
164 /* > LDC >= max(1,N) if NCC > 0;
165 LDC >=1 if NCC = 0. */
166 /* > \endverbatim */
167 /* > */
168 /* > \param[out] WORK */
169 /* > \verbatim */
170 /* > WORK is DOUBLE PRECISION array, dimension (4*N) */
171 /* > \endverbatim */
172 /* > */
173 /* > \param[out] INFO */
174 /* > \verbatim */
175 /* > INFO is INTEGER */
176 /* > = 0: successful exit */
177 /* > < 0: If INFO = -i, the i-th argument had an illegal value */
178 /* > > 0: */
179 /* > if NCVT = NRU = NCC = 0, */
180 /* > = 1, a split was marked by a positive value in E */
181 /* > = 2, current block of Z not diagonalized after 30*N */
182 /* > iterations (in inner while loop) */
183 /* > = 3, termination criterion of outer while loop not met */
184 /* > (program created more than N unreduced blocks) */
185 /* > else NCVT = NRU = NCC = 0, */
186 /* > the algorithm did not converge;
187 D and E contain the */
188 /* > elements of a bidiagonal matrix which is orthogonally */
189 /* > similar to the input matrix B;
190 if INFO = i, i */
191 /* > elements of E have not converged to zero. */
192 /* > \endverbatim */
193 /* > \par Internal Parameters: */
194 /* ========================= */
195 /* > */
196 /* > \verbatim */
197 /* > TOLMUL DOUBLE PRECISION, default = max(10,min(100,EPS**(-1/8))) */
198 /* > TOLMUL controls the convergence criterion of the QR loop. */
199 /* > If it is positive, TOLMUL*EPS is the desired relative */
200 /* > precision in the computed singular values. */
201 /* > If it is negative, f2c_abs(TOLMUL*EPS*sigma_max) is the */
202 /* > desired absolute accuracy in the computed singular */
203 /* > values (corresponds to relative accuracy */
204 /* > f2c_abs(TOLMUL*EPS) in the largest singular value. */
205 /* > f2c_abs(TOLMUL) should be between 1 and 1/EPS, and preferably */
206 /* > between 10 (for fast convergence) and .1/EPS */
207 /* > (for there to be some accuracy in the results). */
208 /* > Default is to lose at either one eighth or 2 of the */
209 /* > available decimal digits in each computed singular value */
210 /* > (whichever is smaller). */
211 /* > */
212 /* > MAXITR INTEGER, default = 6 */
213 /* > MAXITR controls the maximum number of passes of the */
214 /* > algorithm through its inner loop. The algorithms stops */
215 /* > (and so fails to converge) if the number of passes */
216 /* > through the inner loop exceeds MAXITR*N**2. */
217 /* > \endverbatim */
218 /* Authors: */
219 /* ======== */
220 /* > \author Univ. of Tennessee */
221 /* > \author Univ. of California Berkeley */
222 /* > \author Univ. of Colorado Denver */
223 /* > \author NAG Ltd. */
224 /* > \date November 2011 */
225 /* > \ingroup auxOTHERcomputational */
226 /* ===================================================================== */
227 /* Subroutine */
dbdsqr_(char * uplo,integer * n,integer * ncvt,integer * nru,integer * ncc,doublereal * d__,doublereal * e,doublereal * vt,integer * ldvt,doublereal * u,integer * ldu,doublereal * c__,integer * ldc,doublereal * work,integer * info)228 int dbdsqr_(char *uplo, integer *n, integer *ncvt, integer * nru, integer *ncc, doublereal *d__, doublereal *e, doublereal *vt, integer *ldvt, doublereal *u, integer *ldu, doublereal *c__, integer * ldc, doublereal *work, integer *info)
229 {
230 /* System generated locals */
231 integer c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
232 doublereal d__1, d__2, d__3, d__4;
233 /* Builtin functions */
234 double pow_dd(doublereal *, doublereal *), sqrt(doublereal), d_sign( doublereal *, doublereal *);
235 /* Local variables */
236 doublereal f, g, h__;
237 integer i__, j, m;
238 doublereal r__, cs;
239 integer ll;
240 doublereal sn, mu;
241 integer nm1, nm12, nm13, lll;
242 doublereal eps, sll, tol, abse;
243 integer idir;
244 doublereal abss;
245 integer oldm;
246 doublereal cosl;
247 integer isub, iter;
248 doublereal unfl, sinl, cosr, smin, smax, sinr;
249 extern /* Subroutine */
250 int drot_(integer *, doublereal *, integer *, doublereal *, integer *, doublereal *, doublereal *), dlas2_( doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), dscal_(integer *, doublereal *, doublereal *, integer *);
251 extern logical lsame_(char *, char *);
252 doublereal oldcs;
253 extern /* Subroutine */
254 int dlasr_(char *, char *, char *, integer *, integer *, doublereal *, doublereal *, doublereal *, integer *);
255 integer oldll;
256 doublereal shift, sigmn, oldsn;
257 extern /* Subroutine */
258 int dswap_(integer *, doublereal *, integer *, doublereal *, integer *);
259 integer maxit;
260 doublereal sminl, sigmx;
261 logical lower;
262 extern /* Subroutine */
263 int dlasq1_(integer *, doublereal *, doublereal *, doublereal *, integer *), dlasv2_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *, doublereal *);
264 extern doublereal dlamch_(char *);
265 extern /* Subroutine */
266 int dlartg_(doublereal *, doublereal *, doublereal *, doublereal *, doublereal *), xerbla_(char *, integer *);
267 doublereal sminoa, thresh;
268 logical rotate;
269 doublereal tolmul;
270 /* -- LAPACK computational routine (version 3.4.0) -- */
271 /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */
272 /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */
273 /* November 2011 */
274 /* .. Scalar Arguments .. */
275 /* .. */
276 /* .. Array Arguments .. */
277 /* .. */
278 /* ===================================================================== */
279 /* .. Parameters .. */
280 /* .. */
281 /* .. Local Scalars .. */
282 /* .. */
283 /* .. External Functions .. */
284 /* .. */
285 /* .. External Subroutines .. */
286 /* .. */
287 /* .. Intrinsic Functions .. */
288 /* .. */
289 /* .. Executable Statements .. */
290 /* Test the input parameters. */
291 /* Parameter adjustments */
292 --d__;
293 --e;
294 vt_dim1 = *ldvt;
295 vt_offset = 1 + vt_dim1;
296 vt -= vt_offset;
297 u_dim1 = *ldu;
298 u_offset = 1 + u_dim1;
299 u -= u_offset;
300 c_dim1 = *ldc;
301 c_offset = 1 + c_dim1;
302 c__ -= c_offset;
303 --work;
304 /* Function Body */
305 *info = 0;
306 lower = lsame_(uplo, "L");
307 if (! lsame_(uplo, "U") && ! lower)
308 {
309 *info = -1;
310 }
311 else if (*n < 0)
312 {
313 *info = -2;
314 }
315 else if (*ncvt < 0)
316 {
317 *info = -3;
318 }
319 else if (*nru < 0)
320 {
321 *info = -4;
322 }
323 else if (*ncc < 0)
324 {
325 *info = -5;
326 }
327 else if (*ncvt == 0 && *ldvt < 1 || *ncvt > 0 && *ldvt < max(1,*n))
328 {
329 *info = -9;
330 }
331 else if (*ldu < max(1,*nru))
332 {
333 *info = -11;
334 }
335 else if (*ncc == 0 && *ldc < 1 || *ncc > 0 && *ldc < max(1,*n))
336 {
337 *info = -13;
338 }
339 if (*info != 0)
340 {
341 i__1 = -(*info);
342 xerbla_("DBDSQR", &i__1);
343 return 0;
344 }
345 if (*n == 0)
346 {
347 return 0;
348 }
349 if (*n == 1)
350 {
351 goto L160;
352 }
353 /* ROTATE is true if any singular vectors desired, false otherwise */
354 rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;
355 /* If no singular vectors desired, use qd algorithm */
356 if (! rotate)
357 {
358 dlasq1_(n, &d__[1], &e[1], &work[1], info);
359 /* If INFO equals 2, dqds didn't finish, try to finish */
360 if (*info != 2)
361 {
362 return 0;
363 }
364 *info = 0;
365 }
366 nm1 = *n - 1;
367 nm12 = nm1 + nm1;
368 nm13 = nm12 + nm1;
369 idir = 0;
370 /* Get machine constants */
371 eps = dlamch_("Epsilon");
372 unfl = dlamch_("Safe minimum");
373 /* If matrix lower bidiagonal, rotate to be upper bidiagonal */
374 /* by applying Givens rotations on the left */
375 if (lower)
376 {
377 i__1 = *n - 1;
378 for (i__ = 1;
379 i__ <= i__1;
380 ++i__)
381 {
382 dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
383 d__[i__] = r__;
384 e[i__] = sn * d__[i__ + 1];
385 d__[i__ + 1] = cs * d__[i__ + 1];
386 work[i__] = cs;
387 work[nm1 + i__] = sn;
388 /* L10: */
389 }
390 /* Update singular vectors if desired */
391 if (*nru > 0)
392 {
393 dlasr_("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset], ldu);
394 }
395 if (*ncc > 0)
396 {
397 dlasr_("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset], ldc);
398 }
399 }
400 /* Compute singular values to relative accuracy TOL */
401 /* (By setting TOL to be negative, algorithm will compute */
402 /* singular values to absolute accuracy ABS(TOL)*norm(input matrix)) */
403 /* Computing MAX */
404 /* Computing MIN */
405 d__3 = 100.;
406 d__4 = pow_dd(&eps, &c_b15); // , expr subst
407 d__1 = 10.;
408 d__2 = min(d__3,d__4); // , expr subst
409 tolmul = max(d__1,d__2);
410 tol = tolmul * eps;
411 /* Compute approximate maximum, minimum singular values */
412 smax = 0.;
413 i__1 = *n;
414 for (i__ = 1;
415 i__ <= i__1;
416 ++i__)
417 {
418 /* Computing MAX */
419 d__2 = smax;
420 d__3 = (d__1 = d__[i__], f2c_abs(d__1)); // , expr subst
421 smax = max(d__2,d__3);
422 /* L20: */
423 }
424 i__1 = *n - 1;
425 for (i__ = 1;
426 i__ <= i__1;
427 ++i__)
428 {
429 /* Computing MAX */
430 d__2 = smax;
431 d__3 = (d__1 = e[i__], f2c_abs(d__1)); // , expr subst
432 smax = max(d__2,d__3);
433 /* L30: */
434 }
435 sminl = 0.;
436 if (tol >= 0.)
437 {
438 /* Relative accuracy desired */
439 sminoa = f2c_abs(d__[1]);
440 if (sminoa == 0.)
441 {
442 goto L50;
443 }
444 mu = sminoa;
445 i__1 = *n;
446 for (i__ = 2;
447 i__ <= i__1;
448 ++i__)
449 {
450 mu = (d__2 = d__[i__], f2c_abs(d__2)) * (mu / (mu + (d__1 = e[i__ - 1] , f2c_abs(d__1))));
451 sminoa = min(sminoa,mu);
452 if (sminoa == 0.)
453 {
454 goto L50;
455 }
456 /* L40: */
457 }
458 L50:
459 sminoa /= sqrt((doublereal) (*n));
460 /* Computing MAX */
461 d__1 = tol * sminoa;
462 d__2 = *n * 6 * *n * unfl; // , expr subst
463 thresh = max(d__1,d__2);
464 }
465 else
466 {
467 /* Absolute accuracy desired */
468 /* Computing MAX */
469 d__1 = f2c_abs(tol) * smax;
470 d__2 = *n * 6 * *n * unfl; // , expr subst
471 thresh = max(d__1,d__2);
472 }
473 /* Prepare for main iteration loop for the singular values */
474 /* (MAXIT is the maximum number of passes through the inner */
475 /* loop permitted before nonconvergence signalled.) */
476 maxit = *n * 6 * *n;
477 iter = 0;
478 oldll = -1;
479 oldm = -1;
480 /* M points to last element of unconverged part of matrix */
481 m = *n;
482 /* Begin main iteration loop */
483 L60: /* Check for convergence or exceeding iteration count */
484 if (m <= 1)
485 {
486 goto L160;
487 }
488 if (iter > maxit)
489 {
490 goto L200;
491 }
492 /* Find diagonal block of matrix to work on */
493 if (tol < 0. && (d__1 = d__[m], f2c_abs(d__1)) <= thresh)
494 {
495 d__[m] = 0.;
496 }
497 smax = (d__1 = d__[m], f2c_abs(d__1));
498 smin = smax;
499 i__1 = m - 1;
500 for (lll = 1;
501 lll <= i__1;
502 ++lll)
503 {
504 ll = m - lll;
505 abss = (d__1 = d__[ll], f2c_abs(d__1));
506 abse = (d__1 = e[ll], f2c_abs(d__1));
507 if (tol < 0. && abss <= thresh)
508 {
509 d__[ll] = 0.;
510 }
511 if (abse <= thresh)
512 {
513 goto L80;
514 }
515 smin = min(smin,abss);
516 /* Computing MAX */
517 d__1 = max(smax,abss);
518 smax = max(d__1,abse);
519 /* L70: */
520 }
521 ll = 0;
522 goto L90;
523 L80:
524 e[ll] = 0.;
525 /* Matrix splits since E(LL) = 0 */
526 if (ll == m - 1)
527 {
528 /* Convergence of bottom singular value, return to top of loop */
529 --m;
530 goto L60;
531 }
532 L90:
533 ++ll;
534 /* E(LL) through E(M-1) are nonzero, E(LL-1) is zero */
535 if (ll == m - 1)
536 {
537 /* 2 by 2 block, handle separately */
538 dlasv2_(&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr, &sinl, &cosl);
539 d__[m - 1] = sigmx;
540 e[m - 1] = 0.;
541 d__[m] = sigmn;
542 /* Compute singular vectors, if desired */
543 if (*ncvt > 0)
544 {
545 drot_(ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, & cosr, &sinr);
546 }
547 if (*nru > 0)
548 {
549 drot_(nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], & c__1, &cosl, &sinl);
550 }
551 if (*ncc > 0)
552 {
553 drot_(ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, & cosl, &sinl);
554 }
555 m += -2;
556 goto L60;
557 }
558 /* If working on new submatrix, choose shift direction */
559 /* (from larger end diagonal element towards smaller) */
560 if (ll > oldm || m < oldll)
561 {
562 if ((d__1 = d__[ll], f2c_abs(d__1)) >= (d__2 = d__[m], f2c_abs(d__2)))
563 {
564 /* Chase bulge from top (big end) to bottom (small end) */
565 idir = 1;
566 }
567 else
568 {
569 /* Chase bulge from bottom (big end) to top (small end) */
570 idir = 2;
571 }
572 }
573 /* Apply convergence tests */
574 if (idir == 1)
575 {
576 /* Run convergence test in forward direction */
577 /* First apply standard test to bottom of matrix */
578 if ((d__2 = e[m - 1], f2c_abs(d__2)) <= f2c_abs(tol) * (d__1 = d__[m], f2c_abs( d__1)) || tol < 0. && (d__3 = e[m - 1], f2c_abs(d__3)) <= thresh)
579 {
580 e[m - 1] = 0.;
581 goto L60;
582 }
583 if (tol >= 0.)
584 {
585 /* If relative accuracy desired, */
586 /* apply convergence criterion forward */
587 mu = (d__1 = d__[ll], f2c_abs(d__1));
588 sminl = mu;
589 i__1 = m - 1;
590 for (lll = ll;
591 lll <= i__1;
592 ++lll)
593 {
594 if ((d__1 = e[lll], f2c_abs(d__1)) <= tol * mu)
595 {
596 e[lll] = 0.;
597 goto L60;
598 }
599 mu = (d__2 = d__[lll + 1], f2c_abs(d__2)) * (mu / (mu + (d__1 = e[ lll], f2c_abs(d__1))));
600 sminl = min(sminl,mu);
601 /* L100: */
602 }
603 }
604 }
605 else
606 {
607 /* Run convergence test in backward direction */
608 /* First apply standard test to top of matrix */
609 if ((d__2 = e[ll], f2c_abs(d__2)) <= f2c_abs(tol) * (d__1 = d__[ll], f2c_abs(d__1) ) || tol < 0. && (d__3 = e[ll], f2c_abs(d__3)) <= thresh)
610 {
611 e[ll] = 0.;
612 goto L60;
613 }
614 if (tol >= 0.)
615 {
616 /* If relative accuracy desired, */
617 /* apply convergence criterion backward */
618 mu = (d__1 = d__[m], f2c_abs(d__1));
619 sminl = mu;
620 i__1 = ll;
621 for (lll = m - 1;
622 lll >= i__1;
623 --lll)
624 {
625 if ((d__1 = e[lll], f2c_abs(d__1)) <= tol * mu)
626 {
627 e[lll] = 0.;
628 goto L60;
629 }
630 mu = (d__2 = d__[lll], f2c_abs(d__2)) * (mu / (mu + (d__1 = e[lll] , f2c_abs(d__1))));
631 sminl = min(sminl,mu);
632 /* L110: */
633 }
634 }
635 }
636 oldll = ll;
637 oldm = m;
638 /* Compute shift. First, test if shifting would ruin relative */
639 /* accuracy, and if so set the shift to zero. */
640 /* Computing MAX */
641 d__1 = eps;
642 d__2 = tol * .01; // , expr subst
643 if (tol >= 0. && *n * tol * (sminl / smax) <= max(d__1,d__2))
644 {
645 /* Use a zero shift to avoid loss of relative accuracy */
646 shift = 0.;
647 }
648 else
649 {
650 /* Compute the shift from 2-by-2 block at end of matrix */
651 if (idir == 1)
652 {
653 sll = (d__1 = d__[ll], f2c_abs(d__1));
654 dlas2_(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
655 }
656 else
657 {
658 sll = (d__1 = d__[m], f2c_abs(d__1));
659 dlas2_(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
660 }
661 /* Test if shift negligible, and if so set to zero */
662 if (sll > 0.)
663 {
664 /* Computing 2nd power */
665 d__1 = shift / sll;
666 if (d__1 * d__1 < eps)
667 {
668 shift = 0.;
669 }
670 }
671 }
672 /* Increment iteration count */
673 iter = iter + m - ll;
674 /* If SHIFT = 0, do simplified QR iteration */
675 if (shift == 0.)
676 {
677 if (idir == 1)
678 {
679 /* Chase bulge from top to bottom */
680 /* Save cosines and sines for later singular vector updates */
681 cs = 1.;
682 oldcs = 1.;
683 i__1 = m - 1;
684 for (i__ = ll;
685 i__ <= i__1;
686 ++i__)
687 {
688 d__1 = d__[i__] * cs;
689 dlartg_(&d__1, &e[i__], &cs, &sn, &r__);
690 if (i__ > ll)
691 {
692 e[i__ - 1] = oldsn * r__;
693 }
694 d__1 = oldcs * r__;
695 d__2 = d__[i__ + 1] * sn;
696 dlartg_(&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
697 work[i__ - ll + 1] = cs;
698 work[i__ - ll + 1 + nm1] = sn;
699 work[i__ - ll + 1 + nm12] = oldcs;
700 work[i__ - ll + 1 + nm13] = oldsn;
701 /* L120: */
702 }
703 h__ = d__[m] * cs;
704 d__[m] = h__ * oldcs;
705 e[m - 1] = h__ * oldsn;
706 /* Update singular vectors */
707 if (*ncvt > 0)
708 {
709 i__1 = m - ll + 1;
710 dlasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[ ll + vt_dim1], ldvt);
711 }
712 if (*nru > 0)
713 {
714 i__1 = m - ll + 1;
715 dlasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u[ll * u_dim1 + 1], ldu);
716 }
717 if (*ncc > 0)
718 {
719 i__1 = m - ll + 1;
720 dlasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c__[ll + c_dim1], ldc);
721 }
722 /* Test convergence */
723 if ((d__1 = e[m - 1], f2c_abs(d__1)) <= thresh)
724 {
725 e[m - 1] = 0.;
726 }
727 }
728 else
729 {
730 /* Chase bulge from bottom to top */
731 /* Save cosines and sines for later singular vector updates */
732 cs = 1.;
733 oldcs = 1.;
734 i__1 = ll + 1;
735 for (i__ = m;
736 i__ >= i__1;
737 --i__)
738 {
739 d__1 = d__[i__] * cs;
740 dlartg_(&d__1, &e[i__ - 1], &cs, &sn, &r__);
741 if (i__ < m)
742 {
743 e[i__] = oldsn * r__;
744 }
745 d__1 = oldcs * r__;
746 d__2 = d__[i__ - 1] * sn;
747 dlartg_(&d__1, &d__2, &oldcs, &oldsn, &d__[i__]);
748 work[i__ - ll] = cs;
749 work[i__ - ll + nm1] = -sn;
750 work[i__ - ll + nm12] = oldcs;
751 work[i__ - ll + nm13] = -oldsn;
752 /* L130: */
753 }
754 h__ = d__[ll] * cs;
755 d__[ll] = h__ * oldcs;
756 e[ll] = h__ * oldsn;
757 /* Update singular vectors */
758 if (*ncvt > 0)
759 {
760 i__1 = m - ll + 1;
761 dlasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[ nm13 + 1], &vt[ll + vt_dim1], ldvt);
762 }
763 if (*nru > 0)
764 {
765 i__1 = m - ll + 1;
766 dlasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll * u_dim1 + 1], ldu);
767 }
768 if (*ncc > 0)
769 {
770 i__1 = m - ll + 1;
771 dlasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[ ll + c_dim1], ldc);
772 }
773 /* Test convergence */
774 if ((d__1 = e[ll], f2c_abs(d__1)) <= thresh)
775 {
776 e[ll] = 0.;
777 }
778 }
779 }
780 else
781 {
782 /* Use nonzero shift */
783 if (idir == 1)
784 {
785 /* Chase bulge from top to bottom */
786 /* Save cosines and sines for later singular vector updates */
787 f = ((d__1 = d__[ll], f2c_abs(d__1)) - shift) * (d_sign(&c_b49, &d__[ ll]) + shift / d__[ll]);
788 g = e[ll];
789 i__1 = m - 1;
790 for (i__ = ll;
791 i__ <= i__1;
792 ++i__)
793 {
794 dlartg_(&f, &g, &cosr, &sinr, &r__);
795 if (i__ > ll)
796 {
797 e[i__ - 1] = r__;
798 }
799 f = cosr * d__[i__] + sinr * e[i__];
800 e[i__] = cosr * e[i__] - sinr * d__[i__];
801 g = sinr * d__[i__ + 1];
802 d__[i__ + 1] = cosr * d__[i__ + 1];
803 dlartg_(&f, &g, &cosl, &sinl, &r__);
804 d__[i__] = r__;
805 f = cosl * e[i__] + sinl * d__[i__ + 1];
806 d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
807 if (i__ < m - 1)
808 {
809 g = sinl * e[i__ + 1];
810 e[i__ + 1] = cosl * e[i__ + 1];
811 }
812 work[i__ - ll + 1] = cosr;
813 work[i__ - ll + 1 + nm1] = sinr;
814 work[i__ - ll + 1 + nm12] = cosl;
815 work[i__ - ll + 1 + nm13] = sinl;
816 /* L140: */
817 }
818 e[m - 1] = f;
819 /* Update singular vectors */
820 if (*ncvt > 0)
821 {
822 i__1 = m - ll + 1;
823 dlasr_("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[ ll + vt_dim1], ldvt);
824 }
825 if (*nru > 0)
826 {
827 i__1 = m - ll + 1;
828 dlasr_("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u[ll * u_dim1 + 1], ldu);
829 }
830 if (*ncc > 0)
831 {
832 i__1 = m - ll + 1;
833 dlasr_("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c__[ll + c_dim1], ldc);
834 }
835 /* Test convergence */
836 if ((d__1 = e[m - 1], f2c_abs(d__1)) <= thresh)
837 {
838 e[m - 1] = 0.;
839 }
840 }
841 else
842 {
843 /* Chase bulge from bottom to top */
844 /* Save cosines and sines for later singular vector updates */
845 f = ((d__1 = d__[m], f2c_abs(d__1)) - shift) * (d_sign(&c_b49, &d__[m] ) + shift / d__[m]);
846 g = e[m - 1];
847 i__1 = ll + 1;
848 for (i__ = m;
849 i__ >= i__1;
850 --i__)
851 {
852 dlartg_(&f, &g, &cosr, &sinr, &r__);
853 if (i__ < m)
854 {
855 e[i__] = r__;
856 }
857 f = cosr * d__[i__] + sinr * e[i__ - 1];
858 e[i__ - 1] = cosr * e[i__ - 1] - sinr * d__[i__];
859 g = sinr * d__[i__ - 1];
860 d__[i__ - 1] = cosr * d__[i__ - 1];
861 dlartg_(&f, &g, &cosl, &sinl, &r__);
862 d__[i__] = r__;
863 f = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
864 d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
865 if (i__ > ll + 1)
866 {
867 g = sinl * e[i__ - 2];
868 e[i__ - 2] = cosl * e[i__ - 2];
869 }
870 work[i__ - ll] = cosr;
871 work[i__ - ll + nm1] = -sinr;
872 work[i__ - ll + nm12] = cosl;
873 work[i__ - ll + nm13] = -sinl;
874 /* L150: */
875 }
876 e[ll] = f;
877 /* Test convergence */
878 if ((d__1 = e[ll], f2c_abs(d__1)) <= thresh)
879 {
880 e[ll] = 0.;
881 }
882 /* Update singular vectors if desired */
883 if (*ncvt > 0)
884 {
885 i__1 = m - ll + 1;
886 dlasr_("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[ nm13 + 1], &vt[ll + vt_dim1], ldvt);
887 }
888 if (*nru > 0)
889 {
890 i__1 = m - ll + 1;
891 dlasr_("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll * u_dim1 + 1], ldu);
892 }
893 if (*ncc > 0)
894 {
895 i__1 = m - ll + 1;
896 dlasr_("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[ ll + c_dim1], ldc);
897 }
898 }
899 }
900 /* QR iteration finished, go back and check convergence */
901 goto L60;
902 /* All singular values converged, so make them positive */
903 L160:
904 i__1 = *n;
905 for (i__ = 1;
906 i__ <= i__1;
907 ++i__)
908 {
909 if (d__[i__] < 0.)
910 {
911 d__[i__] = -d__[i__];
912 /* Change sign of singular vectors, if desired */
913 if (*ncvt > 0)
914 {
915 dscal_(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
916 }
917 }
918 /* L170: */
919 }
920 /* Sort the singular values into decreasing order (insertion sort on */
921 /* singular values, but only one transposition per singular vector) */
922 i__1 = *n - 1;
923 for (i__ = 1;
924 i__ <= i__1;
925 ++i__)
926 {
927 /* Scan for smallest D(I) */
928 isub = 1;
929 smin = d__[1];
930 i__2 = *n + 1 - i__;
931 for (j = 2;
932 j <= i__2;
933 ++j)
934 {
935 if (d__[j] <= smin)
936 {
937 isub = j;
938 smin = d__[j];
939 }
940 /* L180: */
941 }
942 if (isub != *n + 1 - i__)
943 {
944 /* Swap singular values and vectors */
945 d__[isub] = d__[*n + 1 - i__];
946 d__[*n + 1 - i__] = smin;
947 if (*ncvt > 0)
948 {
949 dswap_(ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ + vt_dim1], ldvt);
950 }
951 if (*nru > 0)
952 {
953 dswap_(nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) * u_dim1 + 1], &c__1);
954 }
955 if (*ncc > 0)
956 {
957 dswap_(ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ + c_dim1], ldc);
958 }
959 }
960 /* L190: */
961 }
962 goto L220;
963 /* Maximum number of iterations exceeded, failure to converge */
964 L200:
965 *info = 0;
966 i__1 = *n - 1;
967 for (i__ = 1;
968 i__ <= i__1;
969 ++i__)
970 {
971 if (e[i__] != 0.)
972 {
973 ++(*info);
974 }
975 /* L210: */
976 }
977 L220:
978 return 0;
979 /* End of DBDSQR */
980 }
981 /* dbdsqr_ */
982