1 /* ./src_f77/sbdsdc.f -- translated by f2c (version 20030320).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include <punc/vf2c.h>
7
8 /* Table of constant values */
9
10 static integer c__9 = 9;
11 static integer c__0 = 0;
12 static real c_b15 = 1.f;
13 static integer c__1 = 1;
14 static real c_b29 = 0.f;
15
sbdsdc_(char * uplo,char * compq,integer * n,real * d__,real * e,real * u,integer * ldu,real * vt,integer * ldvt,real * q,integer * iq,real * work,integer * iwork,integer * info,ftnlen uplo_len,ftnlen compq_len)16 /* Subroutine */ int sbdsdc_(char *uplo, char *compq, integer *n, real *d__,
17 real *e, real *u, integer *ldu, real *vt, integer *ldvt, real *q,
18 integer *iq, real *work, integer *iwork, integer *info, ftnlen
19 uplo_len, ftnlen compq_len)
20 {
21 /* System generated locals */
22 integer u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
23 real r__1;
24
25 /* Builtin functions */
26 double r_sign(real *, real *), log(doublereal);
27
28 /* Local variables */
29 static integer i__, j, k;
30 static real p, r__;
31 static integer z__, ic, ii, kk;
32 static real cs;
33 static integer is, iu;
34 static real sn;
35 static integer nm1;
36 static real eps;
37 static integer ivt, difl, difr, ierr, perm, mlvl, sqre;
38 extern logical lsame_(char *, char *, ftnlen, ftnlen);
39 static integer poles;
40 extern /* Subroutine */ int slasr_(char *, char *, char *, integer *,
41 integer *, real *, real *, real *, integer *, ftnlen, ftnlen,
42 ftnlen);
43 static integer iuplo, nsize, start;
44 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
45 integer *), sswap_(integer *, real *, integer *, real *, integer *
46 ), slasd0_(integer *, integer *, real *, real *, real *, integer *
47 , real *, integer *, integer *, integer *, real *, integer *);
48 extern doublereal slamch_(char *, ftnlen);
49 extern /* Subroutine */ int slasda_(integer *, integer *, integer *,
50 integer *, real *, real *, real *, integer *, real *, integer *,
51 real *, real *, real *, real *, integer *, integer *, integer *,
52 integer *, real *, real *, real *, real *, integer *, integer *),
53 xerbla_(char *, integer *, ftnlen);
54 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
55 integer *, integer *, ftnlen, ftnlen);
56 extern /* Subroutine */ int slascl_(char *, integer *, integer *, real *,
57 real *, integer *, integer *, real *, integer *, integer *,
58 ftnlen);
59 static integer givcol;
60 extern /* Subroutine */ int slasdq_(char *, integer *, integer *, integer
61 *, integer *, integer *, real *, real *, real *, integer *, real *
62 , integer *, real *, integer *, real *, integer *, ftnlen);
63 static integer icompq;
64 extern /* Subroutine */ int slaset_(char *, integer *, integer *, real *,
65 real *, real *, integer *, ftnlen), slartg_(real *, real *, real *
66 , real *, real *);
67 static real orgnrm;
68 static integer givnum;
69 extern doublereal slanst_(char *, integer *, real *, real *, ftnlen);
70 static integer givptr, qstart, smlsiz, wstart, smlszp;
71
72
73 /* -- LAPACK routine (version 3.0) -- */
74 /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
75 /* Courant Institute, Argonne National Lab, and Rice University */
76 /* December 1, 1999 */
77
78 /* .. Scalar Arguments .. */
79 /* .. */
80 /* .. Array Arguments .. */
81 /* .. */
82
83 /* Purpose */
84 /* ======= */
85
86 /* SBDSDC computes the singular value decomposition (SVD) of a real */
87 /* N-by-N (upper or lower) bidiagonal matrix B: B = U * S * VT, */
88 /* using a divide and conquer method, where S is a diagonal matrix */
89 /* with non-negative diagonal elements (the singular values of B), and */
90 /* U and VT are orthogonal matrices of left and right singular vectors, */
91 /* respectively. SBDSDC can be used to compute all singular values, */
92 /* and optionally, singular vectors or singular vectors in compact form. */
93
94 /* This code makes very mild assumptions about floating point */
95 /* arithmetic. It will work on machines with a guard digit in */
96 /* add/subtract, or on those binary machines without guard digits */
97 /* which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. */
98 /* It could conceivably fail on hexadecimal or decimal machines */
99 /* without guard digits, but we know of none. See SLASD3 for details. */
100
101 /* The code currently call SLASDQ if singular values only are desired. */
102 /* However, it can be slightly modified to compute singular values */
103 /* using the divide and conquer method. */
104
105 /* Arguments */
106 /* ========= */
107
108 /* UPLO (input) CHARACTER*1 */
109 /* = 'U': B is upper bidiagonal. */
110 /* = 'L': B is lower bidiagonal. */
111
112 /* COMPQ (input) CHARACTER*1 */
113 /* Specifies whether singular vectors are to be computed */
114 /* as follows: */
115 /* = 'N': Compute singular values only; */
116 /* = 'P': Compute singular values and compute singular */
117 /* vectors in compact form; */
118 /* = 'I': Compute singular values and singular vectors. */
119
120 /* N (input) INTEGER */
121 /* The order of the matrix B. N >= 0. */
122
123 /* D (input/output) REAL array, dimension (N) */
124 /* On entry, the n diagonal elements of the bidiagonal matrix B. */
125 /* On exit, if INFO=0, the singular values of B. */
126
127 /* E (input/output) REAL array, dimension (N) */
128 /* On entry, the elements of E contain the offdiagonal */
129 /* elements of the bidiagonal matrix whose SVD is desired. */
130 /* On exit, E has been destroyed. */
131
132 /* U (output) REAL array, dimension (LDU,N) */
133 /* If COMPQ = 'I', then: */
134 /* On exit, if INFO = 0, U contains the left singular vectors */
135 /* of the bidiagonal matrix. */
136 /* For other values of COMPQ, U is not referenced. */
137
138 /* LDU (input) INTEGER */
139 /* The leading dimension of the array U. LDU >= 1. */
140 /* If singular vectors are desired, then LDU >= max( 1, N ). */
141
142 /* VT (output) REAL array, dimension (LDVT,N) */
143 /* If COMPQ = 'I', then: */
144 /* On exit, if INFO = 0, VT' contains the right singular */
145 /* vectors of the bidiagonal matrix. */
146 /* For other values of COMPQ, VT is not referenced. */
147
148 /* LDVT (input) INTEGER */
149 /* The leading dimension of the array VT. LDVT >= 1. */
150 /* If singular vectors are desired, then LDVT >= max( 1, N ). */
151
152 /* Q (output) REAL array, dimension (LDQ) */
153 /* If COMPQ = 'P', then: */
154 /* On exit, if INFO = 0, Q and IQ contain the left */
155 /* and right singular vectors in a compact form, */
156 /* requiring O(N log N) space instead of 2*N**2. */
157 /* In particular, Q contains all the REAL data in */
158 /* LDQ >= N*(11 + 2*SMLSIZ + 8*INT(LOG_2(N/(SMLSIZ+1)))) */
159 /* words of memory, where SMLSIZ is returned by ILAENV and */
160 /* is equal to the maximum size of the subproblems at the */
161 /* bottom of the computation tree (usually about 25). */
162 /* For other values of COMPQ, Q is not referenced. */
163
164 /* IQ (output) INTEGER array, dimension (LDIQ) */
165 /* If COMPQ = 'P', then: */
166 /* On exit, if INFO = 0, Q and IQ contain the left */
167 /* and right singular vectors in a compact form, */
168 /* requiring O(N log N) space instead of 2*N**2. */
169 /* In particular, IQ contains all INTEGER data in */
170 /* LDIQ >= N*(3 + 3*INT(LOG_2(N/(SMLSIZ+1)))) */
171 /* words of memory, where SMLSIZ is returned by ILAENV and */
172 /* is equal to the maximum size of the subproblems at the */
173 /* bottom of the computation tree (usually about 25). */
174 /* For other values of COMPQ, IQ is not referenced. */
175
176 /* WORK (workspace) REAL array, dimension (LWORK) */
177 /* If COMPQ = 'N' then LWORK >= (4 * N). */
178 /* If COMPQ = 'P' then LWORK >= (6 * N). */
179 /* If COMPQ = 'I' then LWORK >= (3 * N**2 + 4 * N). */
180
181 /* IWORK (workspace) INTEGER array, dimension (8*N) */
182
183 /* INFO (output) INTEGER */
184 /* = 0: successful exit. */
185 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
186 /* > 0: The algorithm failed to compute an singular value. */
187 /* The update process of divide and conquer failed. */
188
189 /* Further Details */
190 /* =============== */
191
192 /* Based on contributions by */
193 /* Ming Gu and Huan Ren, Computer Science Division, University of */
194 /* California at Berkeley, USA */
195
196 /* ===================================================================== */
197
198 /* .. Parameters .. */
199 /* .. */
200 /* .. Local Scalars .. */
201 /* .. */
202 /* .. External Functions .. */
203 /* .. */
204 /* .. External Subroutines .. */
205 /* .. */
206 /* .. Intrinsic Functions .. */
207 /* .. */
208 /* .. Executable Statements .. */
209
210 /* Test the input parameters. */
211
212 /* Parameter adjustments */
213 --d__;
214 --e;
215 u_dim1 = *ldu;
216 u_offset = 1 + u_dim1;
217 u -= u_offset;
218 vt_dim1 = *ldvt;
219 vt_offset = 1 + vt_dim1;
220 vt -= vt_offset;
221 --q;
222 --iq;
223 --work;
224 --iwork;
225
226 /* Function Body */
227 *info = 0;
228
229 iuplo = 0;
230 if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
231 iuplo = 1;
232 }
233 if (lsame_(uplo, "L", (ftnlen)1, (ftnlen)1)) {
234 iuplo = 2;
235 }
236 if (lsame_(compq, "N", (ftnlen)1, (ftnlen)1)) {
237 icompq = 0;
238 } else if (lsame_(compq, "P", (ftnlen)1, (ftnlen)1)) {
239 icompq = 1;
240 } else if (lsame_(compq, "I", (ftnlen)1, (ftnlen)1)) {
241 icompq = 2;
242 } else {
243 icompq = -1;
244 }
245 if (iuplo == 0) {
246 *info = -1;
247 } else if (icompq < 0) {
248 *info = -2;
249 } else if (*n < 0) {
250 *info = -3;
251 } else if (*ldu < 1 || icompq == 2 && *ldu < *n) {
252 *info = -7;
253 } else if (*ldvt < 1 || icompq == 2 && *ldvt < *n) {
254 *info = -9;
255 }
256 if (*info != 0) {
257 i__1 = -(*info);
258 xerbla_("SBDSDC", &i__1, (ftnlen)6);
259 return 0;
260 }
261
262 /* Quick return if possible */
263
264 if (*n == 0) {
265 return 0;
266 }
267 smlsiz = ilaenv_(&c__9, "SBDSDC", " ", &c__0, &c__0, &c__0, &c__0, (
268 ftnlen)6, (ftnlen)1);
269 if (*n == 1) {
270 if (icompq == 1) {
271 q[1] = r_sign(&c_b15, &d__[1]);
272 q[smlsiz * *n + 1] = 1.f;
273 } else if (icompq == 2) {
274 u[u_dim1 + 1] = r_sign(&c_b15, &d__[1]);
275 vt[vt_dim1 + 1] = 1.f;
276 }
277 d__[1] = dabs(d__[1]);
278 return 0;
279 }
280 nm1 = *n - 1;
281
282 /* If matrix lower bidiagonal, rotate to be upper bidiagonal */
283 /* by applying Givens rotations on the left */
284
285 wstart = 1;
286 qstart = 3;
287 if (icompq == 1) {
288 scopy_(n, &d__[1], &c__1, &q[1], &c__1);
289 i__1 = *n - 1;
290 scopy_(&i__1, &e[1], &c__1, &q[*n + 1], &c__1);
291 }
292 if (iuplo == 2) {
293 qstart = 5;
294 wstart = (*n << 1) - 1;
295 i__1 = *n - 1;
296 for (i__ = 1; i__ <= i__1; ++i__) {
297 slartg_(&d__[i__], &e[i__], &cs, &sn, &r__);
298 d__[i__] = r__;
299 e[i__] = sn * d__[i__ + 1];
300 d__[i__ + 1] = cs * d__[i__ + 1];
301 if (icompq == 1) {
302 q[i__ + (*n << 1)] = cs;
303 q[i__ + *n * 3] = sn;
304 } else if (icompq == 2) {
305 work[i__] = cs;
306 work[nm1 + i__] = -sn;
307 }
308 /* L10: */
309 }
310 }
311
312 /* If ICOMPQ = 0, use SLASDQ to compute the singular values. */
313
314 if (icompq == 0) {
315 slasdq_("U", &c__0, n, &c__0, &c__0, &c__0, &d__[1], &e[1], &vt[
316 vt_offset], ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[
317 wstart], info, (ftnlen)1);
318 goto L40;
319 }
320
321 /* If N is smaller than the minimum divide size SMLSIZ, then solve */
322 /* the problem with another solver. */
323
324 if (*n <= smlsiz) {
325 if (icompq == 2) {
326 slaset_("A", n, n, &c_b29, &c_b15, &u[u_offset], ldu, (ftnlen)1);
327 slaset_("A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt, (ftnlen)
328 1);
329 slasdq_("U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &vt[vt_offset]
330 , ldvt, &u[u_offset], ldu, &u[u_offset], ldu, &work[
331 wstart], info, (ftnlen)1);
332 } else if (icompq == 1) {
333 iu = 1;
334 ivt = iu + *n;
335 slaset_("A", n, n, &c_b29, &c_b15, &q[iu + (qstart - 1) * *n], n,
336 (ftnlen)1);
337 slaset_("A", n, n, &c_b29, &c_b15, &q[ivt + (qstart - 1) * *n], n,
338 (ftnlen)1);
339 slasdq_("U", &c__0, n, n, n, &c__0, &d__[1], &e[1], &q[ivt + (
340 qstart - 1) * *n], n, &q[iu + (qstart - 1) * *n], n, &q[
341 iu + (qstart - 1) * *n], n, &work[wstart], info, (ftnlen)
342 1);
343 }
344 goto L40;
345 }
346
347 if (icompq == 2) {
348 slaset_("A", n, n, &c_b29, &c_b15, &u[u_offset], ldu, (ftnlen)1);
349 slaset_("A", n, n, &c_b29, &c_b15, &vt[vt_offset], ldvt, (ftnlen)1);
350 }
351
352 /* Scale. */
353
354 orgnrm = slanst_("M", n, &d__[1], &e[1], (ftnlen)1);
355 if (orgnrm == 0.f) {
356 return 0;
357 }
358 slascl_("G", &c__0, &c__0, &orgnrm, &c_b15, n, &c__1, &d__[1], n, &ierr, (
359 ftnlen)1);
360 slascl_("G", &c__0, &c__0, &orgnrm, &c_b15, &nm1, &c__1, &e[1], &nm1, &
361 ierr, (ftnlen)1);
362
363 eps = slamch_("Epsilon", (ftnlen)7);
364
365 mlvl = (integer) (log((real) (*n) / (real) (smlsiz + 1)) / log(2.f)) + 1;
366 smlszp = smlsiz + 1;
367
368 if (icompq == 1) {
369 iu = 1;
370 ivt = smlsiz + 1;
371 difl = ivt + smlszp;
372 difr = difl + mlvl;
373 z__ = difr + (mlvl << 1);
374 ic = z__ + mlvl;
375 is = ic + 1;
376 poles = is + 1;
377 givnum = poles + (mlvl << 1);
378
379 k = 1;
380 givptr = 2;
381 perm = 3;
382 givcol = perm + mlvl;
383 }
384
385 i__1 = *n;
386 for (i__ = 1; i__ <= i__1; ++i__) {
387 if ((r__1 = d__[i__], dabs(r__1)) < eps) {
388 d__[i__] = r_sign(&eps, &d__[i__]);
389 }
390 /* L20: */
391 }
392
393 start = 1;
394 sqre = 0;
395
396 i__1 = nm1;
397 for (i__ = 1; i__ <= i__1; ++i__) {
398 if ((r__1 = e[i__], dabs(r__1)) < eps || i__ == nm1) {
399
400 /* Subproblem found. First determine its size and then */
401 /* apply divide and conquer on it. */
402
403 if (i__ < nm1) {
404
405 /* A subproblem with E(I) small for I < NM1. */
406
407 nsize = i__ - start + 1;
408 } else if ((r__1 = e[i__], dabs(r__1)) >= eps) {
409
410 /* A subproblem with E(NM1) not too small but I = NM1. */
411
412 nsize = *n - start + 1;
413 } else {
414
415 /* A subproblem with E(NM1) small. This implies an */
416 /* 1-by-1 subproblem at D(N). Solve this 1-by-1 problem */
417 /* first. */
418
419 nsize = i__ - start + 1;
420 if (icompq == 2) {
421 u[*n + *n * u_dim1] = r_sign(&c_b15, &d__[*n]);
422 vt[*n + *n * vt_dim1] = 1.f;
423 } else if (icompq == 1) {
424 q[*n + (qstart - 1) * *n] = r_sign(&c_b15, &d__[*n]);
425 q[*n + (smlsiz + qstart - 1) * *n] = 1.f;
426 }
427 d__[*n] = (r__1 = d__[*n], dabs(r__1));
428 }
429 if (icompq == 2) {
430 slasd0_(&nsize, &sqre, &d__[start], &e[start], &u[start +
431 start * u_dim1], ldu, &vt[start + start * vt_dim1],
432 ldvt, &smlsiz, &iwork[1], &work[wstart], info);
433 } else {
434 slasda_(&icompq, &smlsiz, &nsize, &sqre, &d__[start], &e[
435 start], &q[start + (iu + qstart - 2) * *n], n, &q[
436 start + (ivt + qstart - 2) * *n], &iq[start + k * *n],
437 &q[start + (difl + qstart - 2) * *n], &q[start + (
438 difr + qstart - 2) * *n], &q[start + (z__ + qstart -
439 2) * *n], &q[start + (poles + qstart - 2) * *n], &iq[
440 start + givptr * *n], &iq[start + givcol * *n], n, &
441 iq[start + perm * *n], &q[start + (givnum + qstart -
442 2) * *n], &q[start + (ic + qstart - 2) * *n], &q[
443 start + (is + qstart - 2) * *n], &work[wstart], &
444 iwork[1], info);
445 if (*info != 0) {
446 return 0;
447 }
448 }
449 start = i__ + 1;
450 }
451 /* L30: */
452 }
453
454 /* Unscale */
455
456 slascl_("G", &c__0, &c__0, &c_b15, &orgnrm, n, &c__1, &d__[1], n, &ierr, (
457 ftnlen)1);
458 L40:
459
460 /* Use Selection Sort to minimize swaps of singular vectors */
461
462 i__1 = *n;
463 for (ii = 2; ii <= i__1; ++ii) {
464 i__ = ii - 1;
465 kk = i__;
466 p = d__[i__];
467 i__2 = *n;
468 for (j = ii; j <= i__2; ++j) {
469 if (d__[j] > p) {
470 kk = j;
471 p = d__[j];
472 }
473 /* L50: */
474 }
475 if (kk != i__) {
476 d__[kk] = d__[i__];
477 d__[i__] = p;
478 if (icompq == 1) {
479 iq[i__] = kk;
480 } else if (icompq == 2) {
481 sswap_(n, &u[i__ * u_dim1 + 1], &c__1, &u[kk * u_dim1 + 1], &
482 c__1);
483 sswap_(n, &vt[i__ + vt_dim1], ldvt, &vt[kk + vt_dim1], ldvt);
484 }
485 } else if (icompq == 1) {
486 iq[i__] = i__;
487 }
488 /* L60: */
489 }
490
491 /* If ICOMPQ = 1, use IQ(N,1) as the indicator for UPLO */
492
493 if (icompq == 1) {
494 if (iuplo == 1) {
495 iq[*n] = 1;
496 } else {
497 iq[*n] = 0;
498 }
499 }
500
501 /* If B is lower bidiagonal, update U by those Givens rotations */
502 /* which rotated B to be upper bidiagonal */
503
504 if (iuplo == 2 && icompq == 2) {
505 slasr_("L", "V", "B", n, n, &work[1], &work[*n], &u[u_offset], ldu, (
506 ftnlen)1, (ftnlen)1, (ftnlen)1);
507 }
508
509 return 0;
510
511 /* End of SBDSDC */
512
513 } /* sbdsdc_ */
514
515