1 /* ./src_f77/slasd7.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__1 = 1;
11
slasd7_(integer * icompq,integer * nl,integer * nr,integer * sqre,integer * k,real * d__,real * z__,real * zw,real * vf,real * vfw,real * vl,real * vlw,real * alpha,real * beta,real * dsigma,integer * idx,integer * idxp,integer * idxq,integer * perm,integer * givptr,integer * givcol,integer * ldgcol,real * givnum,integer * ldgnum,real * c__,real * s,integer * info)12 /* Subroutine */ int slasd7_(integer *icompq, integer *nl, integer *nr,
13 integer *sqre, integer *k, real *d__, real *z__, real *zw, real *vf,
14 real *vfw, real *vl, real *vlw, real *alpha, real *beta, real *dsigma,
15 integer *idx, integer *idxp, integer *idxq, integer *perm, integer *
16 givptr, integer *givcol, integer *ldgcol, real *givnum, integer *
17 ldgnum, real *c__, real *s, integer *info)
18 {
19 /* System generated locals */
20 integer givcol_dim1, givcol_offset, givnum_dim1, givnum_offset, i__1;
21 real r__1, r__2;
22
23 /* Local variables */
24 static integer i__, j, m, n, k2;
25 static real z1;
26 static integer jp;
27 static real eps, tau, tol;
28 static integer nlp1, nlp2, idxi, idxj;
29 extern /* Subroutine */ int srot_(integer *, real *, integer *, real *,
30 integer *, real *, real *);
31 static integer idxjp, jprev;
32 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
33 integer *);
34 extern doublereal slapy2_(real *, real *), slamch_(char *, ftnlen);
35 extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), slamrg_(
36 integer *, integer *, real *, integer *, integer *, integer *);
37 static real hlftol;
38
39
40 /* -- LAPACK auxiliary routine (version 3.0) -- */
41 /* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab, */
42 /* Courant Institute, NAG Ltd., and Rice University */
43 /* June 30, 1999 */
44
45 /* .. Scalar Arguments .. */
46 /* .. */
47 /* .. Array Arguments .. */
48 /* .. */
49
50 /* Purpose */
51 /* ======= */
52
53 /* SLASD7 merges the two sets of singular values together into a single */
54 /* sorted set. Then it tries to deflate the size of the problem. There */
55 /* are two ways in which deflation can occur: when two or more singular */
56 /* values are close together or if there is a tiny entry in the Z */
57 /* vector. For each such occurrence the order of the related */
58 /* secular equation problem is reduced by one. */
59
60 /* SLASD7 is called from SLASD6. */
61
62 /* Arguments */
63 /* ========= */
64
65 /* ICOMPQ (input) INTEGER */
66 /* Specifies whether singular vectors are to be computed */
67 /* in compact form, as follows: */
68 /* = 0: Compute singular values only. */
69 /* = 1: Compute singular vectors of upper */
70 /* bidiagonal matrix in compact form. */
71
72 /* NL (input) INTEGER */
73 /* The row dimension of the upper block. NL >= 1. */
74
75 /* NR (input) INTEGER */
76 /* The row dimension of the lower block. NR >= 1. */
77
78 /* SQRE (input) INTEGER */
79 /* = 0: the lower block is an NR-by-NR square matrix. */
80 /* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
81
82 /* The bidiagonal matrix has */
83 /* N = NL + NR + 1 rows and */
84 /* M = N + SQRE >= N columns. */
85
86 /* K (output) INTEGER */
87 /* Contains the dimension of the non-deflated matrix, this is */
88 /* the order of the related secular equation. 1 <= K <=N. */
89
90 /* D (input/output) REAL array, dimension ( N ) */
91 /* On entry D contains the singular values of the two submatrices */
92 /* to be combined. On exit D contains the trailing (N-K) updated */
93 /* singular values (those which were deflated) sorted into */
94 /* increasing order. */
95
96 /* Z (output) REAL array, dimension ( M ) */
97 /* On exit Z contains the updating row vector in the secular */
98 /* equation. */
99
100 /* ZW (workspace) REAL array, dimension ( M ) */
101 /* Workspace for Z. */
102
103 /* VF (input/output) REAL array, dimension ( M ) */
104 /* On entry, VF(1:NL+1) contains the first components of all */
105 /* right singular vectors of the upper block; and VF(NL+2:M) */
106 /* contains the first components of all right singular vectors */
107 /* of the lower block. On exit, VF contains the first components */
108 /* of all right singular vectors of the bidiagonal matrix. */
109
110 /* VFW (workspace) REAL array, dimension ( M ) */
111 /* Workspace for VF. */
112
113 /* VL (input/output) REAL array, dimension ( M ) */
114 /* On entry, VL(1:NL+1) contains the last components of all */
115 /* right singular vectors of the upper block; and VL(NL+2:M) */
116 /* contains the last components of all right singular vectors */
117 /* of the lower block. On exit, VL contains the last components */
118 /* of all right singular vectors of the bidiagonal matrix. */
119
120 /* VLW (workspace) REAL array, dimension ( M ) */
121 /* Workspace for VL. */
122
123 /* ALPHA (input) REAL */
124 /* Contains the diagonal element associated with the added row. */
125
126 /* BETA (input) REAL */
127 /* Contains the off-diagonal element associated with the added */
128 /* row. */
129
130 /* DSIGMA (output) REAL array, dimension ( N ) */
131 /* Contains a copy of the diagonal elements (K-1 singular values */
132 /* and one zero) in the secular equation. */
133
134 /* IDX (workspace) INTEGER array, dimension ( N ) */
135 /* This will contain the permutation used to sort the contents of */
136 /* D into ascending order. */
137
138 /* IDXP (workspace) INTEGER array, dimension ( N ) */
139 /* This will contain the permutation used to place deflated */
140 /* values of D at the end of the array. On output IDXP(2:K) */
141 /* points to the nondeflated D-values and IDXP(K+1:N) */
142 /* points to the deflated singular values. */
143
144 /* IDXQ (input) INTEGER array, dimension ( N ) */
145 /* This contains the permutation which separately sorts the two */
146 /* sub-problems in D into ascending order. Note that entries in */
147 /* the first half of this permutation must first be moved one */
148 /* position backward; and entries in the second half */
149 /* must first have NL+1 added to their values. */
150
151 /* PERM (output) INTEGER array, dimension ( N ) */
152 /* The permutations (from deflation and sorting) to be applied */
153 /* to each singular block. Not referenced if ICOMPQ = 0. */
154
155 /* GIVPTR (output) INTEGER */
156 /* The number of Givens rotations which took place in this */
157 /* subproblem. Not referenced if ICOMPQ = 0. */
158
159 /* GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) */
160 /* Each pair of numbers indicates a pair of columns to take place */
161 /* in a Givens rotation. Not referenced if ICOMPQ = 0. */
162
163 /* LDGCOL (input) INTEGER */
164 /* The leading dimension of GIVCOL, must be at least N. */
165
166 /* GIVNUM (output) REAL array, dimension ( LDGNUM, 2 ) */
167 /* Each number indicates the C or S value to be used in the */
168 /* corresponding Givens rotation. Not referenced if ICOMPQ = 0. */
169
170 /* LDGNUM (input) INTEGER */
171 /* The leading dimension of GIVNUM, must be at least N. */
172
173 /* C (output) REAL */
174 /* C contains garbage if SQRE =0 and the C-value of a Givens */
175 /* rotation related to the right null space if SQRE = 1. */
176
177 /* S (output) REAL */
178 /* S contains garbage if SQRE =0 and the S-value of a Givens */
179 /* rotation related to the right null space if SQRE = 1. */
180
181 /* INFO (output) INTEGER */
182 /* = 0: successful exit. */
183 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
184
185 /* Further Details */
186 /* =============== */
187
188 /* Based on contributions by */
189 /* Ming Gu and Huan Ren, Computer Science Division, University of */
190 /* California at Berkeley, USA */
191
192 /* ===================================================================== */
193
194 /* .. Parameters .. */
195 /* .. */
196 /* .. Local Scalars .. */
197
198 /* .. */
199 /* .. External Subroutines .. */
200 /* .. */
201 /* .. External Functions .. */
202 /* .. */
203 /* .. Intrinsic Functions .. */
204 /* .. */
205 /* .. Executable Statements .. */
206
207 /* Test the input parameters. */
208
209 /* Parameter adjustments */
210 --d__;
211 --z__;
212 --zw;
213 --vf;
214 --vfw;
215 --vl;
216 --vlw;
217 --dsigma;
218 --idx;
219 --idxp;
220 --idxq;
221 --perm;
222 givcol_dim1 = *ldgcol;
223 givcol_offset = 1 + givcol_dim1;
224 givcol -= givcol_offset;
225 givnum_dim1 = *ldgnum;
226 givnum_offset = 1 + givnum_dim1;
227 givnum -= givnum_offset;
228
229 /* Function Body */
230 *info = 0;
231 n = *nl + *nr + 1;
232 m = n + *sqre;
233
234 if (*icompq < 0 || *icompq > 1) {
235 *info = -1;
236 } else if (*nl < 1) {
237 *info = -2;
238 } else if (*nr < 1) {
239 *info = -3;
240 } else if (*sqre < 0 || *sqre > 1) {
241 *info = -4;
242 } else if (*ldgcol < n) {
243 *info = -22;
244 } else if (*ldgnum < n) {
245 *info = -24;
246 }
247 if (*info != 0) {
248 i__1 = -(*info);
249 xerbla_("SLASD7", &i__1, (ftnlen)6);
250 return 0;
251 }
252
253 nlp1 = *nl + 1;
254 nlp2 = *nl + 2;
255 if (*icompq == 1) {
256 *givptr = 0;
257 }
258
259 /* Generate the first part of the vector Z and move the singular */
260 /* values in the first part of D one position backward. */
261
262 z1 = *alpha * vl[nlp1];
263 vl[nlp1] = 0.f;
264 tau = vf[nlp1];
265 for (i__ = *nl; i__ >= 1; --i__) {
266 z__[i__ + 1] = *alpha * vl[i__];
267 vl[i__] = 0.f;
268 vf[i__ + 1] = vf[i__];
269 d__[i__ + 1] = d__[i__];
270 idxq[i__ + 1] = idxq[i__] + 1;
271 /* L10: */
272 }
273 vf[1] = tau;
274
275 /* Generate the second part of the vector Z. */
276
277 i__1 = m;
278 for (i__ = nlp2; i__ <= i__1; ++i__) {
279 z__[i__] = *beta * vf[i__];
280 vf[i__] = 0.f;
281 /* L20: */
282 }
283
284 /* Sort the singular values into increasing order */
285
286 i__1 = n;
287 for (i__ = nlp2; i__ <= i__1; ++i__) {
288 idxq[i__] += nlp1;
289 /* L30: */
290 }
291
292 /* DSIGMA, IDXC, IDXC, and ZW are used as storage space. */
293
294 i__1 = n;
295 for (i__ = 2; i__ <= i__1; ++i__) {
296 dsigma[i__] = d__[idxq[i__]];
297 zw[i__] = z__[idxq[i__]];
298 vfw[i__] = vf[idxq[i__]];
299 vlw[i__] = vl[idxq[i__]];
300 /* L40: */
301 }
302
303 slamrg_(nl, nr, &dsigma[2], &c__1, &c__1, &idx[2]);
304
305 i__1 = n;
306 for (i__ = 2; i__ <= i__1; ++i__) {
307 idxi = idx[i__] + 1;
308 d__[i__] = dsigma[idxi];
309 z__[i__] = zw[idxi];
310 vf[i__] = vfw[idxi];
311 vl[i__] = vlw[idxi];
312 /* L50: */
313 }
314
315 /* Calculate the allowable deflation tolerence */
316
317 eps = slamch_("Epsilon", (ftnlen)7);
318 /* Computing MAX */
319 r__1 = dabs(*alpha), r__2 = dabs(*beta);
320 tol = dmax(r__1,r__2);
321 /* Computing MAX */
322 r__2 = (r__1 = d__[n], dabs(r__1));
323 tol = eps * 64.f * dmax(r__2,tol);
324
325 /* There are 2 kinds of deflation -- first a value in the z-vector */
326 /* is small, second two (or more) singular values are very close */
327 /* together (their difference is small). */
328
329 /* If the value in the z-vector is small, we simply permute the */
330 /* array so that the corresponding singular value is moved to the */
331 /* end. */
332
333 /* If two values in the D-vector are close, we perform a two-sided */
334 /* rotation designed to make one of the corresponding z-vector */
335 /* entries zero, and then permute the array so that the deflated */
336 /* singular value is moved to the end. */
337
338 /* If there are multiple singular values then the problem deflates. */
339 /* Here the number of equal singular values are found. As each equal */
340 /* singular value is found, an elementary reflector is computed to */
341 /* rotate the corresponding singular subspace so that the */
342 /* corresponding components of Z are zero in this new basis. */
343
344 *k = 1;
345 k2 = n + 1;
346 i__1 = n;
347 for (j = 2; j <= i__1; ++j) {
348 if ((r__1 = z__[j], dabs(r__1)) <= tol) {
349
350 /* Deflate due to small z component. */
351
352 --k2;
353 idxp[k2] = j;
354 if (j == n) {
355 goto L100;
356 }
357 } else {
358 jprev = j;
359 goto L70;
360 }
361 /* L60: */
362 }
363 L70:
364 j = jprev;
365 L80:
366 ++j;
367 if (j > n) {
368 goto L90;
369 }
370 if ((r__1 = z__[j], dabs(r__1)) <= tol) {
371
372 /* Deflate due to small z component. */
373
374 --k2;
375 idxp[k2] = j;
376 } else {
377
378 /* Check if singular values are close enough to allow deflation. */
379
380 if ((r__1 = d__[j] - d__[jprev], dabs(r__1)) <= tol) {
381
382 /* Deflation is possible. */
383
384 *s = z__[jprev];
385 *c__ = z__[j];
386
387 /* Find sqrt(a**2+b**2) without overflow or */
388 /* destructive underflow. */
389
390 tau = slapy2_(c__, s);
391 z__[j] = tau;
392 z__[jprev] = 0.f;
393 *c__ /= tau;
394 *s = -(*s) / tau;
395
396 /* Record the appropriate Givens rotation */
397
398 if (*icompq == 1) {
399 ++(*givptr);
400 idxjp = idxq[idx[jprev] + 1];
401 idxj = idxq[idx[j] + 1];
402 if (idxjp <= nlp1) {
403 --idxjp;
404 }
405 if (idxj <= nlp1) {
406 --idxj;
407 }
408 givcol[*givptr + (givcol_dim1 << 1)] = idxjp;
409 givcol[*givptr + givcol_dim1] = idxj;
410 givnum[*givptr + (givnum_dim1 << 1)] = *c__;
411 givnum[*givptr + givnum_dim1] = *s;
412 }
413 srot_(&c__1, &vf[jprev], &c__1, &vf[j], &c__1, c__, s);
414 srot_(&c__1, &vl[jprev], &c__1, &vl[j], &c__1, c__, s);
415 --k2;
416 idxp[k2] = jprev;
417 jprev = j;
418 } else {
419 ++(*k);
420 zw[*k] = z__[jprev];
421 dsigma[*k] = d__[jprev];
422 idxp[*k] = jprev;
423 jprev = j;
424 }
425 }
426 goto L80;
427 L90:
428
429 /* Record the last singular value. */
430
431 ++(*k);
432 zw[*k] = z__[jprev];
433 dsigma[*k] = d__[jprev];
434 idxp[*k] = jprev;
435
436 L100:
437
438 /* Sort the singular values into DSIGMA. The singular values which */
439 /* were not deflated go into the first K slots of DSIGMA, except */
440 /* that DSIGMA(1) is treated separately. */
441
442 i__1 = n;
443 for (j = 2; j <= i__1; ++j) {
444 jp = idxp[j];
445 dsigma[j] = d__[jp];
446 vfw[j] = vf[jp];
447 vlw[j] = vl[jp];
448 /* L110: */
449 }
450 if (*icompq == 1) {
451 i__1 = n;
452 for (j = 2; j <= i__1; ++j) {
453 jp = idxp[j];
454 perm[j] = idxq[idx[jp] + 1];
455 if (perm[j] <= nlp1) {
456 --perm[j];
457 }
458 /* L120: */
459 }
460 }
461
462 /* The deflated singular values go back into the last N - K slots of */
463 /* D. */
464
465 i__1 = n - *k;
466 scopy_(&i__1, &dsigma[*k + 1], &c__1, &d__[*k + 1], &c__1);
467
468 /* Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and */
469 /* VL(M). */
470
471 dsigma[1] = 0.f;
472 hlftol = tol / 2.f;
473 if (dabs(dsigma[2]) <= hlftol) {
474 dsigma[2] = hlftol;
475 }
476 if (m > n) {
477 z__[1] = slapy2_(&z1, &z__[m]);
478 if (z__[1] <= tol) {
479 *c__ = 1.f;
480 *s = 0.f;
481 z__[1] = tol;
482 } else {
483 *c__ = z1 / z__[1];
484 *s = -z__[m] / z__[1];
485 }
486 srot_(&c__1, &vf[m], &c__1, &vf[1], &c__1, c__, s);
487 srot_(&c__1, &vl[m], &c__1, &vl[1], &c__1, c__, s);
488 } else {
489 if (dabs(z1) <= tol) {
490 z__[1] = tol;
491 } else {
492 z__[1] = z1;
493 }
494 }
495
496 /* Restore Z, VF, and VL. */
497
498 i__1 = *k - 1;
499 scopy_(&i__1, &zw[2], &c__1, &z__[2], &c__1);
500 i__1 = n - 1;
501 scopy_(&i__1, &vfw[2], &c__1, &vf[2], &c__1);
502 i__1 = n - 1;
503 scopy_(&i__1, &vlw[2], &c__1, &vl[2], &c__1);
504
505 return 0;
506
507 /* End of SLASD7 */
508
509 } /* slasd7_ */
510
511