1 /* dlasd3.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14 #include "blaswrap.h"
15
16 /* Table of constant values */
17
18 static integer c__1 = 1;
19 static integer c__0 = 0;
20 static doublereal c_b13 = 1.;
21 static doublereal c_b26 = 0.;
22
dlasd3_(integer * nl,integer * nr,integer * sqre,integer * k,doublereal * d__,doublereal * q,integer * ldq,doublereal * dsigma,doublereal * u,integer * ldu,doublereal * u2,integer * ldu2,doublereal * vt,integer * ldvt,doublereal * vt2,integer * ldvt2,integer * idxc,integer * ctot,doublereal * z__,integer * info)23 /* Subroutine */ int dlasd3_(integer *nl, integer *nr, integer *sqre, integer
24 *k, doublereal *d__, doublereal *q, integer *ldq, doublereal *dsigma,
25 doublereal *u, integer *ldu, doublereal *u2, integer *ldu2,
26 doublereal *vt, integer *ldvt, doublereal *vt2, integer *ldvt2,
27 integer *idxc, integer *ctot, doublereal *z__, integer *info)
28 {
29 /* System generated locals */
30 integer q_dim1, q_offset, u_dim1, u_offset, u2_dim1, u2_offset, vt_dim1,
31 vt_offset, vt2_dim1, vt2_offset, i__1, i__2;
32 doublereal d__1, d__2;
33
34 /* Builtin functions */
35 double sqrt(doublereal), d_sign(doublereal *, doublereal *);
36
37 /* Local variables */
38 integer i__, j, m, n, jc;
39 doublereal rho;
40 integer nlp1, nlp2, nrp1;
41 doublereal temp;
42 extern doublereal dnrm2_(integer *, doublereal *, integer *);
43 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
44 integer *, doublereal *, doublereal *, integer *, doublereal *,
45 integer *, doublereal *, doublereal *, integer *);
46 integer ctemp;
47 extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *,
48 doublereal *, integer *);
49 integer ktemp;
50 extern doublereal dlamc3_(doublereal *, doublereal *);
51 extern /* Subroutine */ int dlasd4_(integer *, integer *, doublereal *,
52 doublereal *, doublereal *, doublereal *, doublereal *,
53 doublereal *, integer *), dlascl_(char *, integer *, integer *,
54 doublereal *, doublereal *, integer *, integer *, doublereal *,
55 integer *, integer *), dlacpy_(char *, integer *, integer
56 *, doublereal *, integer *, doublereal *, integer *),
57 xerbla_(char *, integer *);
58
59
60 /* -- LAPACK auxiliary routine (version 3.2) -- */
61 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
62 /* November 2006 */
63
64 /* .. Scalar Arguments .. */
65 /* .. */
66 /* .. Array Arguments .. */
67 /* .. */
68
69 /* Purpose */
70 /* ======= */
71
72 /* DLASD3 finds all the square roots of the roots of the secular */
73 /* equation, as defined by the values in D and Z. It makes the */
74 /* appropriate calls to DLASD4 and then updates the singular */
75 /* vectors by matrix multiplication. */
76
77 /* This code makes very mild assumptions about floating point */
78 /* arithmetic. It will work on machines with a guard digit in */
79 /* add/subtract, or on those binary machines without guard digits */
80 /* which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. */
81 /* It could conceivably fail on hexadecimal or decimal machines */
82 /* without guard digits, but we know of none. */
83
84 /* DLASD3 is called from DLASD1. */
85
86 /* Arguments */
87 /* ========= */
88
89 /* NL (input) INTEGER */
90 /* The row dimension of the upper block. NL >= 1. */
91
92 /* NR (input) INTEGER */
93 /* The row dimension of the lower block. NR >= 1. */
94
95 /* SQRE (input) INTEGER */
96 /* = 0: the lower block is an NR-by-NR square matrix. */
97 /* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. */
98
99 /* The bidiagonal matrix has N = NL + NR + 1 rows and */
100 /* M = N + SQRE >= N columns. */
101
102 /* K (input) INTEGER */
103 /* The size of the secular equation, 1 =< K = < N. */
104
105 /* D (output) DOUBLE PRECISION array, dimension(K) */
106 /* On exit the square roots of the roots of the secular equation, */
107 /* in ascending order. */
108
109 /* Q (workspace) DOUBLE PRECISION array, */
110 /* dimension at least (LDQ,K). */
111
112 /* LDQ (input) INTEGER */
113 /* The leading dimension of the array Q. LDQ >= K. */
114
115 /* DSIGMA (input) DOUBLE PRECISION array, dimension(K) */
116 /* The first K elements of this array contain the old roots */
117 /* of the deflated updating problem. These are the poles */
118 /* of the secular equation. */
119
120 /* U (output) DOUBLE PRECISION array, dimension (LDU, N) */
121 /* The last N - K columns of this matrix contain the deflated */
122 /* left singular vectors. */
123
124 /* LDU (input) INTEGER */
125 /* The leading dimension of the array U. LDU >= N. */
126
127 /* U2 (input/output) DOUBLE PRECISION array, dimension (LDU2, N) */
128 /* The first K columns of this matrix contain the non-deflated */
129 /* left singular vectors for the split problem. */
130
131 /* LDU2 (input) INTEGER */
132 /* The leading dimension of the array U2. LDU2 >= N. */
133
134 /* VT (output) DOUBLE PRECISION array, dimension (LDVT, M) */
135 /* The last M - K columns of VT' contain the deflated */
136 /* right singular vectors. */
137
138 /* LDVT (input) INTEGER */
139 /* The leading dimension of the array VT. LDVT >= N. */
140
141 /* VT2 (input/output) DOUBLE PRECISION array, dimension (LDVT2, N) */
142 /* The first K columns of VT2' contain the non-deflated */
143 /* right singular vectors for the split problem. */
144
145 /* LDVT2 (input) INTEGER */
146 /* The leading dimension of the array VT2. LDVT2 >= N. */
147
148 /* IDXC (input) INTEGER array, dimension ( N ) */
149 /* The permutation used to arrange the columns of U (and rows of */
150 /* VT) into three groups: the first group contains non-zero */
151 /* entries only at and above (or before) NL +1; the second */
152 /* contains non-zero entries only at and below (or after) NL+2; */
153 /* and the third is dense. The first column of U and the row of */
154 /* VT are treated separately, however. */
155
156 /* The rows of the singular vectors found by DLASD4 */
157 /* must be likewise permuted before the matrix multiplies can */
158 /* take place. */
159
160 /* CTOT (input) INTEGER array, dimension ( 4 ) */
161 /* A count of the total number of the various types of columns */
162 /* in U (or rows in VT), as described in IDXC. The fourth column */
163 /* type is any column which has been deflated. */
164
165 /* Z (input) DOUBLE PRECISION array, dimension (K) */
166 /* The first K elements of this array contain the components */
167 /* of the deflation-adjusted updating row vector. */
168
169 /* INFO (output) INTEGER */
170 /* = 0: successful exit. */
171 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
172 /* > 0: if INFO = 1, an singular value did not converge */
173
174 /* Further Details */
175 /* =============== */
176
177 /* Based on contributions by */
178 /* Ming Gu and Huan Ren, Computer Science Division, University of */
179 /* California at Berkeley, USA */
180
181 /* ===================================================================== */
182
183 /* .. Parameters .. */
184 /* .. */
185 /* .. Local Scalars .. */
186 /* .. */
187 /* .. External Functions .. */
188 /* .. */
189 /* .. External Subroutines .. */
190 /* .. */
191 /* .. Intrinsic Functions .. */
192 /* .. */
193 /* .. Executable Statements .. */
194
195 /* Test the input parameters. */
196
197 /* Parameter adjustments */
198 --d__;
199 q_dim1 = *ldq;
200 q_offset = 1 + q_dim1;
201 q -= q_offset;
202 --dsigma;
203 u_dim1 = *ldu;
204 u_offset = 1 + u_dim1;
205 u -= u_offset;
206 u2_dim1 = *ldu2;
207 u2_offset = 1 + u2_dim1;
208 u2 -= u2_offset;
209 vt_dim1 = *ldvt;
210 vt_offset = 1 + vt_dim1;
211 vt -= vt_offset;
212 vt2_dim1 = *ldvt2;
213 vt2_offset = 1 + vt2_dim1;
214 vt2 -= vt2_offset;
215 --idxc;
216 --ctot;
217 --z__;
218
219 /* Function Body */
220 *info = 0;
221
222 if (*nl < 1) {
223 *info = -1;
224 } else if (*nr < 1) {
225 *info = -2;
226 } else if (*sqre != 1 && *sqre != 0) {
227 *info = -3;
228 }
229
230 n = *nl + *nr + 1;
231 m = n + *sqre;
232 nlp1 = *nl + 1;
233 nlp2 = *nl + 2;
234
235 if (*k < 1 || *k > n) {
236 *info = -4;
237 } else if (*ldq < *k) {
238 *info = -7;
239 } else if (*ldu < n) {
240 *info = -10;
241 } else if (*ldu2 < n) {
242 *info = -12;
243 } else if (*ldvt < m) {
244 *info = -14;
245 } else if (*ldvt2 < m) {
246 *info = -16;
247 }
248 if (*info != 0) {
249 i__1 = -(*info);
250 xerbla_("DLASD3", &i__1);
251 return 0;
252 }
253
254 /* Quick return if possible */
255
256 if (*k == 1) {
257 d__[1] = abs(z__[1]);
258 dcopy_(&m, &vt2[vt2_dim1 + 1], ldvt2, &vt[vt_dim1 + 1], ldvt);
259 if (z__[1] > 0.) {
260 dcopy_(&n, &u2[u2_dim1 + 1], &c__1, &u[u_dim1 + 1], &c__1);
261 } else {
262 i__1 = n;
263 for (i__ = 1; i__ <= i__1; ++i__) {
264 u[i__ + u_dim1] = -u2[i__ + u2_dim1];
265 /* L10: */
266 }
267 }
268 return 0;
269 }
270
271 /* Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can */
272 /* be computed with high relative accuracy (barring over/underflow). */
273 /* This is a problem on machines without a guard digit in */
274 /* add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). */
275 /* The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), */
276 /* which on any of these machines zeros out the bottommost */
277 /* bit of DSIGMA(I) if it is 1; this makes the subsequent */
278 /* subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation */
279 /* occurs. On binary machines with a guard digit (almost all */
280 /* machines) it does not change DSIGMA(I) at all. On hexadecimal */
281 /* and decimal machines with a guard digit, it slightly */
282 /* changes the bottommost bits of DSIGMA(I). It does not account */
283 /* for hexadecimal or decimal machines without guard digits */
284 /* (we know of none). We use a subroutine call to compute */
285 /* 2*DSIGMA(I) to prevent optimizing compilers from eliminating */
286 /* this code. */
287
288 i__1 = *k;
289 for (i__ = 1; i__ <= i__1; ++i__) {
290 dsigma[i__] = dlamc3_(&dsigma[i__], &dsigma[i__]) - dsigma[i__];
291 /* L20: */
292 }
293
294 /* Keep a copy of Z. */
295
296 dcopy_(k, &z__[1], &c__1, &q[q_offset], &c__1);
297
298 /* Normalize Z. */
299
300 rho = dnrm2_(k, &z__[1], &c__1);
301 dlascl_("G", &c__0, &c__0, &rho, &c_b13, k, &c__1, &z__[1], k, info);
302 rho *= rho;
303
304 /* Find the new singular values. */
305
306 i__1 = *k;
307 for (j = 1; j <= i__1; ++j) {
308 dlasd4_(k, &j, &dsigma[1], &z__[1], &u[j * u_dim1 + 1], &rho, &d__[j],
309 &vt[j * vt_dim1 + 1], info);
310
311 /* If the zero finder fails, the computation is terminated. */
312
313 if (*info != 0) {
314 return 0;
315 }
316 /* L30: */
317 }
318
319 /* Compute updated Z. */
320
321 i__1 = *k;
322 for (i__ = 1; i__ <= i__1; ++i__) {
323 z__[i__] = u[i__ + *k * u_dim1] * vt[i__ + *k * vt_dim1];
324 i__2 = i__ - 1;
325 for (j = 1; j <= i__2; ++j) {
326 z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
327 i__] - dsigma[j]) / (dsigma[i__] + dsigma[j]);
328 /* L40: */
329 }
330 i__2 = *k - 1;
331 for (j = i__; j <= i__2; ++j) {
332 z__[i__] *= u[i__ + j * u_dim1] * vt[i__ + j * vt_dim1] / (dsigma[
333 i__] - dsigma[j + 1]) / (dsigma[i__] + dsigma[j + 1]);
334 /* L50: */
335 }
336 d__2 = sqrt((d__1 = z__[i__], abs(d__1)));
337 z__[i__] = d_sign(&d__2, &q[i__ + q_dim1]);
338 /* L60: */
339 }
340
341 /* Compute left singular vectors of the modified diagonal matrix, */
342 /* and store related information for the right singular vectors. */
343
344 i__1 = *k;
345 for (i__ = 1; i__ <= i__1; ++i__) {
346 vt[i__ * vt_dim1 + 1] = z__[1] / u[i__ * u_dim1 + 1] / vt[i__ *
347 vt_dim1 + 1];
348 u[i__ * u_dim1 + 1] = -1.;
349 i__2 = *k;
350 for (j = 2; j <= i__2; ++j) {
351 vt[j + i__ * vt_dim1] = z__[j] / u[j + i__ * u_dim1] / vt[j + i__
352 * vt_dim1];
353 u[j + i__ * u_dim1] = dsigma[j] * vt[j + i__ * vt_dim1];
354 /* L70: */
355 }
356 temp = dnrm2_(k, &u[i__ * u_dim1 + 1], &c__1);
357 q[i__ * q_dim1 + 1] = u[i__ * u_dim1 + 1] / temp;
358 i__2 = *k;
359 for (j = 2; j <= i__2; ++j) {
360 jc = idxc[j];
361 q[j + i__ * q_dim1] = u[jc + i__ * u_dim1] / temp;
362 /* L80: */
363 }
364 /* L90: */
365 }
366
367 /* Update the left singular vector matrix. */
368
369 if (*k == 2) {
370 dgemm_("N", "N", &n, k, k, &c_b13, &u2[u2_offset], ldu2, &q[q_offset],
371 ldq, &c_b26, &u[u_offset], ldu);
372 goto L100;
373 }
374 if (ctot[1] > 0) {
375 dgemm_("N", "N", nl, k, &ctot[1], &c_b13, &u2[(u2_dim1 << 1) + 1],
376 ldu2, &q[q_dim1 + 2], ldq, &c_b26, &u[u_dim1 + 1], ldu);
377 if (ctot[3] > 0) {
378 ktemp = ctot[1] + 2 + ctot[2];
379 dgemm_("N", "N", nl, k, &ctot[3], &c_b13, &u2[ktemp * u2_dim1 + 1]
380 , ldu2, &q[ktemp + q_dim1], ldq, &c_b13, &u[u_dim1 + 1],
381 ldu);
382 }
383 } else if (ctot[3] > 0) {
384 ktemp = ctot[1] + 2 + ctot[2];
385 dgemm_("N", "N", nl, k, &ctot[3], &c_b13, &u2[ktemp * u2_dim1 + 1],
386 ldu2, &q[ktemp + q_dim1], ldq, &c_b26, &u[u_dim1 + 1], ldu);
387 } else {
388 dlacpy_("F", nl, k, &u2[u2_offset], ldu2, &u[u_offset], ldu);
389 }
390 dcopy_(k, &q[q_dim1 + 1], ldq, &u[nlp1 + u_dim1], ldu);
391 ktemp = ctot[1] + 2;
392 ctemp = ctot[2] + ctot[3];
393 dgemm_("N", "N", nr, k, &ctemp, &c_b13, &u2[nlp2 + ktemp * u2_dim1], ldu2,
394 &q[ktemp + q_dim1], ldq, &c_b26, &u[nlp2 + u_dim1], ldu);
395
396 /* Generate the right singular vectors. */
397
398 L100:
399 i__1 = *k;
400 for (i__ = 1; i__ <= i__1; ++i__) {
401 temp = dnrm2_(k, &vt[i__ * vt_dim1 + 1], &c__1);
402 q[i__ + q_dim1] = vt[i__ * vt_dim1 + 1] / temp;
403 i__2 = *k;
404 for (j = 2; j <= i__2; ++j) {
405 jc = idxc[j];
406 q[i__ + j * q_dim1] = vt[jc + i__ * vt_dim1] / temp;
407 /* L110: */
408 }
409 /* L120: */
410 }
411
412 /* Update the right singular vector matrix. */
413
414 if (*k == 2) {
415 dgemm_("N", "N", k, &m, k, &c_b13, &q[q_offset], ldq, &vt2[vt2_offset]
416 , ldvt2, &c_b26, &vt[vt_offset], ldvt);
417 return 0;
418 }
419 ktemp = ctot[1] + 1;
420 dgemm_("N", "N", k, &nlp1, &ktemp, &c_b13, &q[q_dim1 + 1], ldq, &vt2[
421 vt2_dim1 + 1], ldvt2, &c_b26, &vt[vt_dim1 + 1], ldvt);
422 ktemp = ctot[1] + 2 + ctot[2];
423 if (ktemp <= *ldvt2) {
424 dgemm_("N", "N", k, &nlp1, &ctot[3], &c_b13, &q[ktemp * q_dim1 + 1],
425 ldq, &vt2[ktemp + vt2_dim1], ldvt2, &c_b13, &vt[vt_dim1 + 1],
426 ldvt);
427 }
428
429 ktemp = ctot[1] + 1;
430 nrp1 = *nr + *sqre;
431 if (ktemp > 1) {
432 i__1 = *k;
433 for (i__ = 1; i__ <= i__1; ++i__) {
434 q[i__ + ktemp * q_dim1] = q[i__ + q_dim1];
435 /* L130: */
436 }
437 i__1 = m;
438 for (i__ = nlp2; i__ <= i__1; ++i__) {
439 vt2[ktemp + i__ * vt2_dim1] = vt2[i__ * vt2_dim1 + 1];
440 /* L140: */
441 }
442 }
443 ctemp = ctot[2] + 1 + ctot[3];
444 dgemm_("N", "N", k, &nrp1, &ctemp, &c_b13, &q[ktemp * q_dim1 + 1], ldq, &
445 vt2[ktemp + nlp2 * vt2_dim1], ldvt2, &c_b26, &vt[nlp2 * vt_dim1 +
446 1], ldvt);
447
448 return 0;
449
450 /* End of DLASD3 */
451
452 } /* dlasd3_ */
453