3 /* Table of constant values */
5 static integer c__1 = 1;
6 static integer c__2 = 2;
7 static integer c__0 = 0;
9 /* Subroutine */ int dlasq1_(integer *n, doublereal *d__, doublereal *e,
10 doublereal *work, integer *info)
12 /* System generated locals */
14 doublereal d__1, d__2, d__3;
16 /* Builtin functions */
17 double sqrt(doublereal);
22 extern /* Subroutine */ int dlas2_(doublereal *, doublereal *, doublereal
23 *, doublereal *, doublereal *);
27 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
28 doublereal *, integer *);
30 extern /* Subroutine */ int dlasq2_(integer *, doublereal *, integer *);
31 extern doublereal dlamch_(char *);
32 extern /* Subroutine */ int dlascl_(char *, integer *, integer *,
33 doublereal *, doublereal *, integer *, integer *, doublereal *,
34 integer *, integer *);
36 extern /* Subroutine */ int xerbla_(char *, integer *), dlasrt_(
37 char *, integer *, doublereal *, integer *);
40 /* -- LAPACK routine (version 3.1) -- */
41 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
44 /* .. Scalar Arguments .. */
46 /* .. Array Arguments .. */
52 /* DLASQ1 computes the singular values of a real N-by-N bidiagonal */
53 /* matrix with diagonal D and off-diagonal E. The singular values */
54 /* are computed to high relative accuracy, in the absence of */
55 /* denormalization, underflow and overflow. The algorithm was first */
58 /* "Accurate singular values and differential qd algorithms" by K. V. */
59 /* Fernando and B. N. Parlett, Numer. Math., Vol-67, No. 2, pp. 191-230, */
62 /* and the present implementation is described in "An implementation of */
63 /* the dqds Algorithm (Positive Case)", LAPACK Working Note. */
68 /* N (input) INTEGER */
69 /* The number of rows and columns in the matrix. N >= 0. */
71 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
72 /* On entry, D contains the diagonal elements of the */
73 /* bidiagonal matrix whose SVD is desired. On normal exit, */
74 /* D contains the singular values in decreasing order. */
76 /* E (input/output) DOUBLE PRECISION array, dimension (N) */
77 /* On entry, elements E(1:N-1) contain the off-diagonal elements */
78 /* of the bidiagonal matrix whose SVD is desired. */
79 /* On exit, E is overwritten. */
81 /* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
83 /* INFO (output) INTEGER */
84 /* = 0: successful exit */
85 /* < 0: if INFO = -i, the i-th argument had an illegal value */
86 /* > 0: the algorithm failed */
87 /* = 1, a split was marked by a positive value in E */
88 /* = 2, current block of Z not diagonalized after 30*N */
89 /* iterations (in inner while loop) */
90 /* = 3, termination criterion of outer while loop not met */
91 /* (program created more than N unreduced blocks) */
93 /* ===================================================================== */
95 /* .. Parameters .. */
97 /* .. Local Scalars .. */
99 /* .. External Subroutines .. */
101 /* .. External Functions .. */
103 /* .. Intrinsic Functions .. */
105 /* .. Executable Statements .. */
107 /* Parameter adjustments */
117 xerbla_("DLASQ1", &i__1);
119 } else if (*n == 0) {
121 } else if (*n == 1) {
122 d__[1] = abs(d__[1]);
124 } else if (*n == 2) {
125 dlas2_(&d__[1], &e[1], &d__[2], &sigmn, &sigmx);
131 /* Estimate the largest singular value. */
135 for (i__ = 1; i__ <= i__1; ++i__) {
136 d__[i__] = (d__1 = d__[i__], abs(d__1));
138 d__2 = sigmx, d__3 = (d__1 = e[i__], abs(d__1));
139 sigmx = max(d__2,d__3);
142 d__[*n] = (d__1 = d__[*n], abs(d__1));
144 /* Early return if SIGMX is zero (matrix is already diagonal). */
147 dlasrt_("D", n, &d__[1], &iinfo);
152 for (i__ = 1; i__ <= i__1; ++i__) {
154 d__1 = sigmx, d__2 = d__[i__];
155 sigmx = max(d__1,d__2);
159 /* Copy D and E into WORK (in the Z format) and scale (squaring the */
160 /* input data makes scaling by a power of the radix pointless). */
162 eps = dlamch_("Precision");
163 safmin = dlamch_("Safe minimum");
164 scale = sqrt(eps / safmin);
165 dcopy_(n, &d__[1], &c__1, &work[1], &c__2);
167 dcopy_(&i__1, &e[1], &c__1, &work[2], &c__2);
168 i__1 = (*n << 1) - 1;
169 i__2 = (*n << 1) - 1;
170 dlascl_("G", &c__0, &c__0, &sigmx, &scale, &i__1, &c__1, &work[1], &i__2,
173 /* Compute the q's and e's. */
175 i__1 = (*n << 1) - 1;
176 for (i__ = 1; i__ <= i__1; ++i__) {
177 /* Computing 2nd power */
179 work[i__] = d__1 * d__1;
184 dlasq2_(n, &work[1], info);
188 for (i__ = 1; i__ <= i__1; ++i__) {
189 d__[i__] = sqrt(work[i__]);
192 dlascl_("G", &c__0, &c__0, &scale, &sigmx, n, &c__1, &d__[1], n, &