3 /* Table of constant values */
5 static real c_b7 = 1.f;
6 static real c_b8 = 0.f;
7 static integer c__2 = 2;
9 /* Subroutine */ int slalsa_(integer *icompq, integer *smlsiz, integer *n,
10 integer *nrhs, real *b, integer *ldb, real *bx, integer *ldbx, real *
11 u, integer *ldu, real *vt, integer *k, real *difl, real *difr, real *
12 z__, real *poles, integer *givptr, integer *givcol, integer *ldgcol,
13 integer *perm, real *givnum, real *c__, real *s, real *work, integer *
16 /* System generated locals */
17 integer givcol_dim1, givcol_offset, perm_dim1, perm_offset, b_dim1,
18 b_offset, bx_dim1, bx_offset, difl_dim1, difl_offset, difr_dim1,
19 difr_offset, givnum_dim1, givnum_offset, poles_dim1, poles_offset,
20 u_dim1, u_offset, vt_dim1, vt_offset, z_dim1, z_offset, i__1,
23 /* Builtin functions */
24 integer pow_ii(integer *, integer *);
27 integer i__, j, i1, ic, lf, nd, ll, nl, nr, im1, nlf, nrf, lvl, ndb1,
28 nlp1, lvl2, nrp1, nlvl, sqre, inode, ndiml;
29 extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
30 integer *, real *, real *, integer *, real *, integer *, real *,
33 extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *,
34 integer *), slals0_(integer *, integer *, integer *, integer *,
35 integer *, real *, integer *, real *, integer *, integer *,
36 integer *, integer *, integer *, real *, integer *, real *, real *
37 , real *, real *, integer *, real *, real *, real *, integer *),
38 xerbla_(char *, integer *), slasdt_(integer *, integer *,
39 integer *, integer *, integer *, integer *, integer *);
42 /* -- LAPACK routine (version 3.1) -- */
43 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
46 /* .. Scalar Arguments .. */
48 /* .. Array Arguments .. */
54 /* SLALSA is an itermediate step in solving the least squares problem */
55 /* by computing the SVD of the coefficient matrix in compact form (The */
56 /* singular vectors are computed as products of simple orthorgonal */
59 /* If ICOMPQ = 0, SLALSA applies the inverse of the left singular vector */
60 /* matrix of an upper bidiagonal matrix to the right hand side; and if */
61 /* ICOMPQ = 1, SLALSA applies the right singular vector matrix to the */
62 /* right hand side. The singular vector matrices were generated in */
63 /* compact form by SLALSA. */
69 /* ICOMPQ (input) INTEGER */
70 /* Specifies whether the left or the right singular vector */
71 /* matrix is involved. */
72 /* = 0: Left singular vector matrix */
73 /* = 1: Right singular vector matrix */
75 /* SMLSIZ (input) INTEGER */
76 /* The maximum size of the subproblems at the bottom of the */
77 /* computation tree. */
79 /* N (input) INTEGER */
80 /* The row and column dimensions of the upper bidiagonal matrix. */
82 /* NRHS (input) INTEGER */
83 /* The number of columns of B and BX. NRHS must be at least 1. */
85 /* B (input/output) REAL array, dimension ( LDB, NRHS ) */
86 /* On input, B contains the right hand sides of the least */
87 /* squares problem in rows 1 through M. */
88 /* On output, B contains the solution X in rows 1 through N. */
90 /* LDB (input) INTEGER */
91 /* The leading dimension of B in the calling subprogram. */
92 /* LDB must be at least max(1,MAX( M, N ) ). */
94 /* BX (output) REAL array, dimension ( LDBX, NRHS ) */
95 /* On exit, the result of applying the left or right singular */
96 /* vector matrix to B. */
98 /* LDBX (input) INTEGER */
99 /* The leading dimension of BX. */
101 /* U (input) REAL array, dimension ( LDU, SMLSIZ ). */
102 /* On entry, U contains the left singular vector matrices of all */
103 /* subproblems at the bottom level. */
105 /* LDU (input) INTEGER, LDU = > N. */
106 /* The leading dimension of arrays U, VT, DIFL, DIFR, */
107 /* POLES, GIVNUM, and Z. */
109 /* VT (input) REAL array, dimension ( LDU, SMLSIZ+1 ). */
110 /* On entry, VT' contains the right singular vector matrices of */
111 /* all subproblems at the bottom level. */
113 /* K (input) INTEGER array, dimension ( N ). */
115 /* DIFL (input) REAL array, dimension ( LDU, NLVL ). */
116 /* where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1. */
118 /* DIFR (input) REAL array, dimension ( LDU, 2 * NLVL ). */
119 /* On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record */
120 /* distances between singular values on the I-th level and */
121 /* singular values on the (I -1)-th level, and DIFR(*, 2 * I) */
122 /* record the normalizing factors of the right singular vectors */
123 /* matrices of subproblems on I-th level. */
125 /* Z (input) REAL array, dimension ( LDU, NLVL ). */
126 /* On entry, Z(1, I) contains the components of the deflation- */
127 /* adjusted updating row vector for subproblems on the I-th */
130 /* POLES (input) REAL array, dimension ( LDU, 2 * NLVL ). */
131 /* On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old */
132 /* singular values involved in the secular equations on the I-th */
135 /* GIVPTR (input) INTEGER array, dimension ( N ). */
136 /* On entry, GIVPTR( I ) records the number of Givens */
137 /* rotations performed on the I-th problem on the computation */
140 /* GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 * NLVL ). */
141 /* On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the */
142 /* locations of Givens rotations performed on the I-th level on */
143 /* the computation tree. */
145 /* LDGCOL (input) INTEGER, LDGCOL = > N. */
146 /* The leading dimension of arrays GIVCOL and PERM. */
148 /* PERM (input) INTEGER array, dimension ( LDGCOL, NLVL ). */
149 /* On entry, PERM(*, I) records permutations done on the I-th */
150 /* level of the computation tree. */
152 /* GIVNUM (input) REAL array, dimension ( LDU, 2 * NLVL ). */
153 /* On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S- */
154 /* values of Givens rotations performed on the I-th level on the */
155 /* computation tree. */
157 /* C (input) REAL array, dimension ( N ). */
158 /* On entry, if the I-th subproblem is not square, */
159 /* C( I ) contains the C-value of a Givens rotation related to */
160 /* the right null space of the I-th subproblem. */
162 /* S (input) REAL array, dimension ( N ). */
163 /* On entry, if the I-th subproblem is not square, */
164 /* S( I ) contains the S-value of a Givens rotation related to */
165 /* the right null space of the I-th subproblem. */
167 /* WORK (workspace) REAL array. */
168 /* The dimension must be at least N. */
170 /* IWORK (workspace) INTEGER array. */
171 /* The dimension must be at least 3 * N */
173 /* INFO (output) INTEGER */
174 /* = 0: successful exit. */
175 /* < 0: if INFO = -i, the i-th argument had an illegal value. */
177 /* Further Details */
178 /* =============== */
180 /* Based on contributions by */
181 /* Ming Gu and Ren-Cang Li, Computer Science Division, University of */
182 /* California at Berkeley, USA */
183 /* Osni Marques, LBNL/NERSC, USA */
185 /* ===================================================================== */
187 /* .. Parameters .. */
189 /* .. Local Scalars .. */
191 /* .. External Subroutines .. */
193 /* .. Executable Statements .. */
195 /* Test the input parameters. */
197 /* Parameter adjustments */
199 b_offset = 1 + b_dim1;
202 bx_offset = 1 + bx_dim1;
205 givnum_offset = 1 + givnum_dim1;
206 givnum -= givnum_offset;
208 poles_offset = 1 + poles_dim1;
209 poles -= poles_offset;
211 z_offset = 1 + z_dim1;
214 difr_offset = 1 + difr_dim1;
217 difl_offset = 1 + difl_dim1;
220 vt_offset = 1 + vt_dim1;
223 u_offset = 1 + u_dim1;
228 perm_offset = 1 + perm_dim1;
230 givcol_dim1 = *ldgcol;
231 givcol_offset = 1 + givcol_dim1;
232 givcol -= givcol_offset;
241 if (*icompq < 0 || *icompq > 1) {
243 } else if (*smlsiz < 3) {
245 } else if (*n < *smlsiz) {
247 } else if (*nrhs < 1) {
249 } else if (*ldb < *n) {
251 } else if (*ldbx < *n) {
253 } else if (*ldu < *n) {
255 } else if (*ldgcol < *n) {
260 xerbla_("SLALSA", &i__1);
264 /* Book-keeping and setting up the computation tree. */
270 slasdt_(n, &nlvl, &nd, &iwork[inode], &iwork[ndiml], &iwork[ndimr],
273 /* The following code applies back the left singular vector factors. */
274 /* For applying back the right singular vector factors, go to 50. */
280 /* The nodes on the bottom level of the tree were solved */
281 /* by SLASDQ. The corresponding left and right singular vector */
282 /* matrices are in explicit form. First apply back the left */
283 /* singular vector matrices. */
287 for (i__ = ndb1; i__ <= i__1; ++i__) {
289 /* IC : center row of each node */
290 /* NL : number of rows of left subproblem */
291 /* NR : number of rows of right subproblem */
292 /* NLF: starting row of the left subproblem */
293 /* NRF: starting row of the right subproblem */
296 ic = iwork[inode + i1];
297 nl = iwork[ndiml + i1];
298 nr = iwork[ndimr + i1];
301 sgemm_("T", "N", &nl, nrhs, &nl, &c_b7, &u[nlf + u_dim1], ldu, &b[nlf
302 + b_dim1], ldb, &c_b8, &bx[nlf + bx_dim1], ldbx);
303 sgemm_("T", "N", &nr, nrhs, &nr, &c_b7, &u[nrf + u_dim1], ldu, &b[nrf
304 + b_dim1], ldb, &c_b8, &bx[nrf + bx_dim1], ldbx);
308 /* Next copy the rows of B that correspond to unchanged rows */
309 /* in the bidiagonal matrix to BX. */
312 for (i__ = 1; i__ <= i__1; ++i__) {
313 ic = iwork[inode + i__ - 1];
314 scopy_(nrhs, &b[ic + b_dim1], ldb, &bx[ic + bx_dim1], ldbx);
318 /* Finally go through the left singular vector matrices of all */
319 /* the other subproblems bottom-up on the tree. */
321 j = pow_ii(&c__2, &nlvl);
324 for (lvl = nlvl; lvl >= 1; --lvl) {
325 lvl2 = (lvl << 1) - 1;
327 /* find the first node LF and last node LL on */
328 /* the current level LVL */
335 lf = pow_ii(&c__2, &i__1);
339 for (i__ = lf; i__ <= i__1; ++i__) {
341 ic = iwork[inode + im1];
342 nl = iwork[ndiml + im1];
343 nr = iwork[ndimr + im1];
347 slals0_(icompq, &nl, &nr, &sqre, nrhs, &bx[nlf + bx_dim1], ldbx, &
348 b[nlf + b_dim1], ldb, &perm[nlf + lvl * perm_dim1], &
349 givptr[j], &givcol[nlf + lvl2 * givcol_dim1], ldgcol, &
350 givnum[nlf + lvl2 * givnum_dim1], ldu, &poles[nlf + lvl2 *
351 poles_dim1], &difl[nlf + lvl * difl_dim1], &difr[nlf +
352 lvl2 * difr_dim1], &z__[nlf + lvl * z_dim1], &k[j], &c__[
353 j], &s[j], &work[1], info);
360 /* ICOMPQ = 1: applying back the right singular vector factors. */
364 /* First now go through the right singular vector matrices of all */
365 /* the tree nodes top-down. */
369 for (lvl = 1; lvl <= i__1; ++lvl) {
370 lvl2 = (lvl << 1) - 1;
372 /* Find the first node LF and last node LL on */
373 /* the current level LVL. */
380 lf = pow_ii(&c__2, &i__2);
384 for (i__ = ll; i__ >= i__2; --i__) {
386 ic = iwork[inode + im1];
387 nl = iwork[ndiml + im1];
388 nr = iwork[ndimr + im1];
397 slals0_(icompq, &nl, &nr, &sqre, nrhs, &b[nlf + b_dim1], ldb, &bx[
398 nlf + bx_dim1], ldbx, &perm[nlf + lvl * perm_dim1], &
399 givptr[j], &givcol[nlf + lvl2 * givcol_dim1], ldgcol, &
400 givnum[nlf + lvl2 * givnum_dim1], ldu, &poles[nlf + lvl2 *
401 poles_dim1], &difl[nlf + lvl * difl_dim1], &difr[nlf +
402 lvl2 * difr_dim1], &z__[nlf + lvl * z_dim1], &k[j], &c__[
403 j], &s[j], &work[1], info);
409 /* The nodes on the bottom level of the tree were solved */
410 /* by SLASDQ. The corresponding right singular vector */
411 /* matrices are in explicit form. Apply them back. */
415 for (i__ = ndb1; i__ <= i__1; ++i__) {
417 ic = iwork[inode + i1];
418 nl = iwork[ndiml + i1];
419 nr = iwork[ndimr + i1];
428 sgemm_("T", "N", &nlp1, nrhs, &nlp1, &c_b7, &vt[nlf + vt_dim1], ldu, &
429 b[nlf + b_dim1], ldb, &c_b8, &bx[nlf + bx_dim1], ldbx);
430 sgemm_("T", "N", &nrp1, nrhs, &nrp1, &c_b7, &vt[nrf + vt_dim1], ldu, &
431 b[nrf + b_dim1], ldb, &c_b8, &bx[nrf + bx_dim1], ldbx);