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