3 /* Table of constant values */
5 static integer c__1 = 1;
6 static doublereal c_b6 = 0.;
7 static integer c__0 = 0;
8 static doublereal c_b11 = 1.;
10 /* Subroutine */ int dlalsd_(char *uplo, integer *smlsiz, integer *n, integer
11 *nrhs, doublereal *d__, doublereal *e, doublereal *b, integer *ldb,
12 doublereal *rcond, integer *rank, doublereal *work, integer *iwork,
15 /* System generated locals */
16 integer b_dim1, b_offset, i__1, i__2;
19 /* Builtin functions */
20 double log(doublereal), d_sign(doublereal *, doublereal *);
23 integer c__, i__, j, k;
29 integer st, vt, nm1, st1;
36 extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
37 doublereal *, integer *, doublereal *, doublereal *);
38 integer nlvl, sqre, bxst;
39 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
40 integer *, doublereal *, doublereal *, integer *, doublereal *,
41 integer *, doublereal *, doublereal *, integer *),
42 dcopy_(integer *, doublereal *, integer *, doublereal *, integer
44 integer poles, sizei, nsize, nwork, icmpq1, icmpq2;
45 extern doublereal dlamch_(char *);
46 extern /* Subroutine */ int dlasda_(integer *, integer *, integer *,
47 integer *, doublereal *, doublereal *, doublereal *, integer *,
48 doublereal *, integer *, doublereal *, doublereal *, doublereal *,
49 doublereal *, integer *, integer *, integer *, integer *,
50 doublereal *, doublereal *, doublereal *, doublereal *, integer *,
51 integer *), dlalsa_(integer *, integer *, integer *, integer *,
52 doublereal *, integer *, doublereal *, integer *, doublereal *,
53 integer *, doublereal *, integer *, doublereal *, doublereal *,
54 doublereal *, doublereal *, integer *, integer *, integer *,
55 integer *, doublereal *, doublereal *, doublereal *, doublereal *,
56 integer *, integer *), dlascl_(char *, integer *, integer *,
57 doublereal *, doublereal *, integer *, integer *, doublereal *,
58 integer *, integer *);
59 extern integer idamax_(integer *, doublereal *, integer *);
60 extern /* Subroutine */ int dlasdq_(char *, integer *, integer *, integer
61 *, integer *, integer *, doublereal *, doublereal *, doublereal *,
62 integer *, doublereal *, integer *, doublereal *, integer *,
63 doublereal *, integer *), dlacpy_(char *, integer *,
64 integer *, doublereal *, integer *, doublereal *, integer *), dlartg_(doublereal *, doublereal *, doublereal *,
65 doublereal *, doublereal *), dlaset_(char *, integer *, integer *,
66 doublereal *, doublereal *, doublereal *, integer *),
67 xerbla_(char *, integer *);
69 extern doublereal dlanst_(char *, integer *, doublereal *, doublereal *);
70 extern /* Subroutine */ int dlasrt_(char *, integer *, doublereal *,
73 integer givnum, givptr, smlszp;
76 /* -- LAPACK routine (version 3.1) -- */
77 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
80 /* .. Scalar Arguments .. */
82 /* .. Array Arguments .. */
88 /* DLALSD uses the singular value decomposition of A to solve the least */
89 /* squares problem of finding X to minimize the Euclidean norm of each */
90 /* column of A*X-B, where A is N-by-N upper bidiagonal, and X and B */
91 /* are N-by-NRHS. The solution X overwrites B. */
93 /* The singular values of A smaller than RCOND times the largest */
94 /* singular value are treated as zero in solving the least squares */
95 /* problem; in this case a minimum norm solution is returned. */
96 /* The actual singular values are returned in D in ascending order. */
98 /* This code makes very mild assumptions about floating point */
99 /* arithmetic. It will work on machines with a guard digit in */
100 /* add/subtract, or on those binary machines without guard digits */
101 /* which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. */
102 /* It could conceivably fail on hexadecimal or decimal machines */
103 /* without guard digits, but we know of none. */
108 /* UPLO (input) CHARACTER*1 */
109 /* = 'U': D and E define an upper bidiagonal matrix. */
110 /* = 'L': D and E define a lower bidiagonal matrix. */
112 /* SMLSIZ (input) INTEGER */
113 /* The maximum size of the subproblems at the bottom of the */
114 /* computation tree. */
116 /* N (input) INTEGER */
117 /* The dimension of the bidiagonal matrix. N >= 0. */
119 /* NRHS (input) INTEGER */
120 /* The number of columns of B. NRHS must be at least 1. */
122 /* D (input/output) DOUBLE PRECISION array, dimension (N) */
123 /* On entry D contains the main diagonal of the bidiagonal */
124 /* matrix. On exit, if INFO = 0, D contains its singular values. */
126 /* E (input/output) DOUBLE PRECISION array, dimension (N-1) */
127 /* Contains the super-diagonal entries of the bidiagonal matrix. */
128 /* On exit, E has been destroyed. */
130 /* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
131 /* On input, B contains the right hand sides of the least */
132 /* squares problem. On output, B contains the solution X. */
134 /* LDB (input) INTEGER */
135 /* The leading dimension of B in the calling subprogram. */
136 /* LDB must be at least max(1,N). */
138 /* RCOND (input) DOUBLE PRECISION */
139 /* The singular values of A less than or equal to RCOND times */
140 /* the largest singular value are treated as zero in solving */
141 /* the least squares problem. If RCOND is negative, */
142 /* machine precision is used instead. */
143 /* For example, if diag(S)*X=B were the least squares problem, */
144 /* where diag(S) is a diagonal matrix of singular values, the */
145 /* solution would be X(i) = B(i) / S(i) if S(i) is greater than */
146 /* RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to */
149 /* RANK (output) INTEGER */
150 /* The number of singular values of A greater than RCOND times */
151 /* the largest singular value. */
153 /* WORK (workspace) DOUBLE PRECISION array, dimension at least */
154 /* (9*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2), */
155 /* where NLVL = max(0, INT(log_2 (N/(SMLSIZ+1))) + 1). */
157 /* IWORK (workspace) INTEGER array, dimension at least */
158 /* (3*N*NLVL + 11*N) */
160 /* INFO (output) INTEGER */
161 /* = 0: successful exit. */
162 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
163 /* > 0: The algorithm failed to compute an singular value while */
164 /* working on the submatrix lying in rows and columns */
165 /* INFO/(N+1) through MOD(INFO,N+1). */
167 /* Further Details */
168 /* =============== */
170 /* Based on contributions by */
171 /* Ming Gu and Ren-Cang Li, Computer Science Division, University of */
172 /* California at Berkeley, USA */
173 /* Osni Marques, LBNL/NERSC, USA */
175 /* ===================================================================== */
177 /* .. Parameters .. */
179 /* .. Local Scalars .. */
181 /* .. External Functions .. */
183 /* .. External Subroutines .. */
185 /* .. Intrinsic Functions .. */
187 /* .. Executable Statements .. */
189 /* Test the input parameters. */
191 /* Parameter adjustments */
195 b_offset = 1 + b_dim1;
205 } else if (*nrhs < 1) {
207 } else if (*ldb < 1 || *ldb < *n) {
212 xerbla_("DLALSD", &i__1);
216 eps = dlamch_("Epsilon");
218 /* Set up the tolerance. */
220 if (*rcond <= 0. || *rcond >= 1.) {
228 /* Quick return if possible. */
232 } else if (*n == 1) {
234 dlaset_("A", &c__1, nrhs, &c_b6, &c_b6, &b[b_offset], ldb);
237 dlascl_("G", &c__0, &c__0, &d__[1], &c_b11, &c__1, nrhs, &b[
238 b_offset], ldb, info);
239 d__[1] = abs(d__[1]);
244 /* Rotate the matrix if it is lower bidiagonal. */
246 if (*(unsigned char *)uplo == 'L') {
248 for (i__ = 1; i__ <= i__1; ++i__) {
249 dlartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
251 e[i__] = sn * d__[i__ + 1];
252 d__[i__ + 1] = cs * d__[i__ + 1];
254 drot_(&c__1, &b[i__ + b_dim1], &c__1, &b[i__ + 1 + b_dim1], &
257 work[(i__ << 1) - 1] = cs;
264 for (i__ = 1; i__ <= i__1; ++i__) {
266 for (j = 1; j <= i__2; ++j) {
267 cs = work[(j << 1) - 1];
269 drot_(&c__1, &b[j + i__ * b_dim1], &c__1, &b[j + 1 + i__ *
270 b_dim1], &c__1, &cs, &sn);
281 orgnrm = dlanst_("M", n, &d__[1], &e[1]);
283 dlaset_("A", n, nrhs, &c_b6, &c_b6, &b[b_offset], ldb);
287 dlascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, &c__1, &d__[1], n, info);
288 dlascl_("G", &c__0, &c__0, &orgnrm, &c_b11, &nm1, &c__1, &e[1], &nm1,
291 /* If N is smaller than the minimum divide size SMLSIZ, then solve */
292 /* the problem with another solver. */
296 dlaset_("A", n, n, &c_b6, &c_b11, &work[1], n);
297 dlasdq_("U", &c__0, n, n, &c__0, nrhs, &d__[1], &e[1], &work[1], n, &
298 work[1], n, &b[b_offset], ldb, &work[nwork], info);
302 tol = rcnd * (d__1 = d__[idamax_(n, &d__[1], &c__1)], abs(d__1));
304 for (i__ = 1; i__ <= i__1; ++i__) {
305 if (d__[i__] <= tol) {
306 dlaset_("A", &c__1, nrhs, &c_b6, &c_b6, &b[i__ + b_dim1], ldb);
308 dlascl_("G", &c__0, &c__0, &d__[i__], &c_b11, &c__1, nrhs, &b[
309 i__ + b_dim1], ldb, info);
314 dgemm_("T", "N", n, nrhs, n, &c_b11, &work[1], n, &b[b_offset], ldb, &
315 c_b6, &work[nwork], n);
316 dlacpy_("A", n, nrhs, &work[nwork], n, &b[b_offset], ldb);
320 dlascl_("G", &c__0, &c__0, &c_b11, &orgnrm, n, &c__1, &d__[1], n,
322 dlasrt_("D", n, &d__[1], info);
323 dlascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, nrhs, &b[b_offset],
329 /* Book-keeping and setting up some constants. */
331 nlvl = (integer) (log((doublereal) (*n) / (doublereal) (*smlsiz + 1)) /
334 smlszp = *smlsiz + 1;
337 vt = *smlsiz * *n + 1;
338 difl = vt + smlszp * *n;
339 difr = difl + nlvl * *n;
340 z__ = difr + (nlvl * *n << 1);
341 c__ = z__ + nlvl * *n;
344 givnum = poles + (nlvl << 1) * *n;
345 bx = givnum + (nlvl << 1) * *n;
346 nwork = bx + *n * *nrhs;
352 givcol = perm + nlvl * *n;
353 iwk = givcol + (nlvl * *n << 1);
362 for (i__ = 1; i__ <= i__1; ++i__) {
363 if ((d__1 = d__[i__], abs(d__1)) < eps) {
364 d__[i__] = d_sign(&eps, &d__[i__]);
370 for (i__ = 1; i__ <= i__1; ++i__) {
371 if ((d__1 = e[i__], abs(d__1)) < eps || i__ == nm1) {
375 /* Subproblem found. First determine its size and then */
376 /* apply divide and conquer on it. */
380 /* A subproblem with E(I) small for I < NM1. */
382 nsize = i__ - st + 1;
383 iwork[sizei + nsub - 1] = nsize;
384 } else if ((d__1 = e[i__], abs(d__1)) >= eps) {
386 /* A subproblem with E(NM1) not too small but I = NM1. */
389 iwork[sizei + nsub - 1] = nsize;
392 /* A subproblem with E(NM1) small. This implies an */
393 /* 1-by-1 subproblem at D(N), which is not solved */
396 nsize = i__ - st + 1;
397 iwork[sizei + nsub - 1] = nsize;
400 iwork[sizei + nsub - 1] = 1;
401 dcopy_(nrhs, &b[*n + b_dim1], ldb, &work[bx + nm1], n);
406 /* This is a 1-by-1 subproblem and is not solved */
409 dcopy_(nrhs, &b[st + b_dim1], ldb, &work[bx + st1], n);
410 } else if (nsize <= *smlsiz) {
412 /* This is a small subproblem and is solved by DLASDQ. */
414 dlaset_("A", &nsize, &nsize, &c_b6, &c_b11, &work[vt + st1],
416 dlasdq_("U", &c__0, &nsize, &nsize, &c__0, nrhs, &d__[st], &e[
417 st], &work[vt + st1], n, &work[nwork], n, &b[st +
418 b_dim1], ldb, &work[nwork], info);
422 dlacpy_("A", &nsize, nrhs, &b[st + b_dim1], ldb, &work[bx +
426 /* A large problem. Solve it using divide and conquer. */
428 dlasda_(&icmpq1, smlsiz, &nsize, &sqre, &d__[st], &e[st], &
429 work[u + st1], n, &work[vt + st1], &iwork[k + st1], &
430 work[difl + st1], &work[difr + st1], &work[z__ + st1],
431 &work[poles + st1], &iwork[givptr + st1], &iwork[
432 givcol + st1], n, &iwork[perm + st1], &work[givnum +
433 st1], &work[c__ + st1], &work[s + st1], &work[nwork],
439 dlalsa_(&icmpq2, smlsiz, &nsize, nrhs, &b[st + b_dim1], ldb, &
440 work[bxst], n, &work[u + st1], n, &work[vt + st1], &
441 iwork[k + st1], &work[difl + st1], &work[difr + st1],
442 &work[z__ + st1], &work[poles + st1], &iwork[givptr +
443 st1], &iwork[givcol + st1], n, &iwork[perm + st1], &
444 work[givnum + st1], &work[c__ + st1], &work[s + st1],
445 &work[nwork], &iwork[iwk], info);
455 /* Apply the singular values and treat the tiny ones as zero. */
457 tol = rcnd * (d__1 = d__[idamax_(n, &d__[1], &c__1)], abs(d__1));
460 for (i__ = 1; i__ <= i__1; ++i__) {
462 /* Some of the elements in D can be negative because 1-by-1 */
463 /* subproblems were not solved explicitly. */
465 if ((d__1 = d__[i__], abs(d__1)) <= tol) {
466 dlaset_("A", &c__1, nrhs, &c_b6, &c_b6, &work[bx + i__ - 1], n);
469 dlascl_("G", &c__0, &c__0, &d__[i__], &c_b11, &c__1, nrhs, &work[
470 bx + i__ - 1], n, info);
472 d__[i__] = (d__1 = d__[i__], abs(d__1));
476 /* Now apply back the right singular vectors. */
480 for (i__ = 1; i__ <= i__1; ++i__) {
483 nsize = iwork[sizei + i__ - 1];
486 dcopy_(nrhs, &work[bxst], n, &b[st + b_dim1], ldb);
487 } else if (nsize <= *smlsiz) {
488 dgemm_("T", "N", &nsize, nrhs, &nsize, &c_b11, &work[vt + st1], n,
489 &work[bxst], n, &c_b6, &b[st + b_dim1], ldb);
491 dlalsa_(&icmpq2, smlsiz, &nsize, nrhs, &work[bxst], n, &b[st +
492 b_dim1], ldb, &work[u + st1], n, &work[vt + st1], &iwork[
493 k + st1], &work[difl + st1], &work[difr + st1], &work[z__
494 + st1], &work[poles + st1], &iwork[givptr + st1], &iwork[
495 givcol + st1], n, &iwork[perm + st1], &work[givnum + st1],
496 &work[c__ + st1], &work[s + st1], &work[nwork], &iwork[
505 /* Unscale and sort the singular values. */
507 dlascl_("G", &c__0, &c__0, &c_b11, &orgnrm, n, &c__1, &d__[1], n, info);
508 dlasrt_("D", n, &d__[1], info);
509 dlascl_("G", &c__0, &c__0, &orgnrm, &c_b11, n, nrhs, &b[b_offset], ldb,