LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dlalsa.f
Go to the documentation of this file.
1 *> \brief \b DLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLALSA + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlalsa.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlalsa.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlalsa.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DLALSA( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
22 * LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR,
23 * GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK,
24 * IWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * INTEGER ICOMPQ, INFO, LDB, LDBX, LDGCOL, LDU, N, NRHS,
28 * $ SMLSIZ
29 * ..
30 * .. Array Arguments ..
31 * INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
32 * $ K( * ), PERM( LDGCOL, * )
33 * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), C( * ),
34 * $ DIFL( LDU, * ), DIFR( LDU, * ),
35 * $ GIVNUM( LDU, * ), POLES( LDU, * ), S( * ),
36 * $ U( LDU, * ), VT( LDU, * ), WORK( * ),
37 * $ Z( LDU, * )
38 * ..
39 *
40 *
41 *> \par Purpose:
42 * =============
43 *>
44 *> \verbatim
45 *>
46 *> DLALSA is an itermediate step in solving the least squares problem
47 *> by computing the SVD of the coefficient matrix in compact form (The
48 *> singular vectors are computed as products of simple orthorgonal
49 *> matrices.).
50 *>
51 *> If ICOMPQ = 0, DLALSA applies the inverse of the left singular vector
52 *> matrix of an upper bidiagonal matrix to the right hand side; and if
53 *> ICOMPQ = 1, DLALSA applies the right singular vector matrix to the
54 *> right hand side. The singular vector matrices were generated in
55 *> compact form by DLALSA.
56 *> \endverbatim
57 *
58 * Arguments:
59 * ==========
60 *
61 *> \param[in] ICOMPQ
62 *> \verbatim
63 *> ICOMPQ is INTEGER
64 *> Specifies whether the left or the right singular vector
65 *> matrix is involved.
66 *> = 0: Left singular vector matrix
67 *> = 1: Right singular vector matrix
68 *> \endverbatim
69 *>
70 *> \param[in] SMLSIZ
71 *> \verbatim
72 *> SMLSIZ is INTEGER
73 *> The maximum size of the subproblems at the bottom of the
74 *> computation tree.
75 *> \endverbatim
76 *>
77 *> \param[in] N
78 *> \verbatim
79 *> N is INTEGER
80 *> The row and column dimensions of the upper bidiagonal matrix.
81 *> \endverbatim
82 *>
83 *> \param[in] NRHS
84 *> \verbatim
85 *> NRHS is INTEGER
86 *> The number of columns of B and BX. NRHS must be at least 1.
87 *> \endverbatim
88 *>
89 *> \param[in,out] B
90 *> \verbatim
91 *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS )
92 *> On input, B contains the right hand sides of the least
93 *> squares problem in rows 1 through M.
94 *> On output, B contains the solution X in rows 1 through N.
95 *> \endverbatim
96 *>
97 *> \param[in] LDB
98 *> \verbatim
99 *> LDB is INTEGER
100 *> The leading dimension of B in the calling subprogram.
101 *> LDB must be at least max(1,MAX( M, N ) ).
102 *> \endverbatim
103 *>
104 *> \param[out] BX
105 *> \verbatim
106 *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS )
107 *> On exit, the result of applying the left or right singular
108 *> vector matrix to B.
109 *> \endverbatim
110 *>
111 *> \param[in] LDBX
112 *> \verbatim
113 *> LDBX is INTEGER
114 *> The leading dimension of BX.
115 *> \endverbatim
116 *>
117 *> \param[in] U
118 *> \verbatim
119 *> U is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ ).
120 *> On entry, U contains the left singular vector matrices of all
121 *> subproblems at the bottom level.
122 *> \endverbatim
123 *>
124 *> \param[in] LDU
125 *> \verbatim
126 *> LDU is INTEGER, LDU = > N.
127 *> The leading dimension of arrays U, VT, DIFL, DIFR,
128 *> POLES, GIVNUM, and Z.
129 *> \endverbatim
130 *>
131 *> \param[in] VT
132 *> \verbatim
133 *> VT is DOUBLE PRECISION array, dimension ( LDU, SMLSIZ+1 ).
134 *> On entry, VT**T contains the right singular vector matrices of
135 *> all subproblems at the bottom level.
136 *> \endverbatim
137 *>
138 *> \param[in] K
139 *> \verbatim
140 *> K is INTEGER array, dimension ( N ).
141 *> \endverbatim
142 *>
143 *> \param[in] DIFL
144 *> \verbatim
145 *> DIFL is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
146 *> where NLVL = INT(log_2 (N/(SMLSIZ+1))) + 1.
147 *> \endverbatim
148 *>
149 *> \param[in] DIFR
150 *> \verbatim
151 *> DIFR is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
152 *> On entry, DIFL(*, I) and DIFR(*, 2 * I -1) record
153 *> distances between singular values on the I-th level and
154 *> singular values on the (I -1)-th level, and DIFR(*, 2 * I)
155 *> record the normalizing factors of the right singular vectors
156 *> matrices of subproblems on I-th level.
157 *> \endverbatim
158 *>
159 *> \param[in] Z
160 *> \verbatim
161 *> Z is DOUBLE PRECISION array, dimension ( LDU, NLVL ).
162 *> On entry, Z(1, I) contains the components of the deflation-
163 *> adjusted updating row vector for subproblems on the I-th
164 *> level.
165 *> \endverbatim
166 *>
167 *> \param[in] POLES
168 *> \verbatim
169 *> POLES is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
170 *> On entry, POLES(*, 2 * I -1: 2 * I) contains the new and old
171 *> singular values involved in the secular equations on the I-th
172 *> level.
173 *> \endverbatim
174 *>
175 *> \param[in] GIVPTR
176 *> \verbatim
177 *> GIVPTR is INTEGER array, dimension ( N ).
178 *> On entry, GIVPTR( I ) records the number of Givens
179 *> rotations performed on the I-th problem on the computation
180 *> tree.
181 *> \endverbatim
182 *>
183 *> \param[in] GIVCOL
184 *> \verbatim
185 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 * NLVL ).
186 *> On entry, for each I, GIVCOL(*, 2 * I - 1: 2 * I) records the
187 *> locations of Givens rotations performed on the I-th level on
188 *> the computation tree.
189 *> \endverbatim
190 *>
191 *> \param[in] LDGCOL
192 *> \verbatim
193 *> LDGCOL is INTEGER, LDGCOL = > N.
194 *> The leading dimension of arrays GIVCOL and PERM.
195 *> \endverbatim
196 *>
197 *> \param[in] PERM
198 *> \verbatim
199 *> PERM is INTEGER array, dimension ( LDGCOL, NLVL ).
200 *> On entry, PERM(*, I) records permutations done on the I-th
201 *> level of the computation tree.
202 *> \endverbatim
203 *>
204 *> \param[in] GIVNUM
205 *> \verbatim
206 *> GIVNUM is DOUBLE PRECISION array, dimension ( LDU, 2 * NLVL ).
207 *> On entry, GIVNUM(*, 2 *I -1 : 2 * I) records the C- and S-
208 *> values of Givens rotations performed on the I-th level on the
209 *> computation tree.
210 *> \endverbatim
211 *>
212 *> \param[in] C
213 *> \verbatim
214 *> C is DOUBLE PRECISION array, dimension ( N ).
215 *> On entry, if the I-th subproblem is not square,
216 *> C( I ) contains the C-value of a Givens rotation related to
217 *> the right null space of the I-th subproblem.
218 *> \endverbatim
219 *>
220 *> \param[in] S
221 *> \verbatim
222 *> S is DOUBLE PRECISION array, dimension ( N ).
223 *> On entry, if the I-th subproblem is not square,
224 *> S( I ) contains the S-value of a Givens rotation related to
225 *> the right null space of the I-th subproblem.
226 *> \endverbatim
227 *>
228 *> \param[out] WORK
229 *> \verbatim
230 *> WORK is DOUBLE PRECISION array.
231 *> The dimension must be at least N.
232 *> \endverbatim
233 *>
234 *> \param[out] IWORK
235 *> \verbatim
236 *> IWORK is INTEGER array.
237 *> The dimension must be at least 3 * N
238 *> \endverbatim
239 *>
240 *> \param[out] INFO
241 *> \verbatim
242 *> INFO is INTEGER
243 *> = 0: successful exit.
244 *> < 0: if INFO = -i, the i-th argument had an illegal value.
245 *> \endverbatim
246 *
247 * Authors:
248 * ========
249 *
250 *> \author Univ. of Tennessee
251 *> \author Univ. of California Berkeley
252 *> \author Univ. of Colorado Denver
253 *> \author NAG Ltd.
254 *
255 *> \date September 2012
256 *
257 *> \ingroup doubleOTHERcomputational
258 *
259 *> \par Contributors:
260 * ==================
261 *>
262 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
263 *> California at Berkeley, USA \n
264 *> Osni Marques, LBNL/NERSC, USA \n
265 *
266 * =====================================================================
267  SUBROUTINE dlalsa( ICOMPQ, SMLSIZ, N, NRHS, B, LDB, BX, LDBX, U,
268  $ ldu, vt, k, difl, difr, z, poles, givptr,
269  $ givcol, ldgcol, perm, givnum, c, s, work,
270  $ iwork, info )
271 *
272 * -- LAPACK computational routine (version 3.4.2) --
273 * -- LAPACK is a software package provided by Univ. of Tennessee, --
274 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
275 * September 2012
276 *
277 * .. Scalar Arguments ..
278  INTEGER icompq, info, ldb, ldbx, ldgcol, ldu, n, nrhs,
279  $ smlsiz
280 * ..
281 * .. Array Arguments ..
282  INTEGER givcol( ldgcol, * ), givptr( * ), iwork( * ),
283  $ k( * ), perm( ldgcol, * )
284  DOUBLE PRECISION b( ldb, * ), bx( ldbx, * ), c( * ),
285  $ difl( ldu, * ), difr( ldu, * ),
286  $ givnum( ldu, * ), poles( ldu, * ), s( * ),
287  $ u( ldu, * ), vt( ldu, * ), work( * ),
288  $ z( ldu, * )
289 * ..
290 *
291 * =====================================================================
292 *
293 * .. Parameters ..
294  DOUBLE PRECISION zero, one
295  parameter( zero = 0.0d0, one = 1.0d0 )
296 * ..
297 * .. Local Scalars ..
298  INTEGER i, i1, ic, im1, inode, j, lf, ll, lvl, lvl2,
299  $ nd, ndb1, ndiml, ndimr, nl, nlf, nlp1, nlvl,
300  $ nr, nrf, nrp1, sqre
301 * ..
302 * .. External Subroutines ..
303  EXTERNAL dcopy, dgemm, dlals0, dlasdt, xerbla
304 * ..
305 * .. Executable Statements ..
306 *
307 * Test the input parameters.
308 *
309  info = 0
310 *
311  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
312  info = -1
313  ELSE IF( smlsiz.LT.3 ) THEN
314  info = -2
315  ELSE IF( n.LT.smlsiz ) THEN
316  info = -3
317  ELSE IF( nrhs.LT.1 ) THEN
318  info = -4
319  ELSE IF( ldb.LT.n ) THEN
320  info = -6
321  ELSE IF( ldbx.LT.n ) THEN
322  info = -8
323  ELSE IF( ldu.LT.n ) THEN
324  info = -10
325  ELSE IF( ldgcol.LT.n ) THEN
326  info = -19
327  END IF
328  IF( info.NE.0 ) THEN
329  CALL xerbla( 'DLALSA', -info )
330  return
331  END IF
332 *
333 * Book-keeping and setting up the computation tree.
334 *
335  inode = 1
336  ndiml = inode + n
337  ndimr = ndiml + n
338 *
339  CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
340  $ iwork( ndimr ), smlsiz )
341 *
342 * The following code applies back the left singular vector factors.
343 * For applying back the right singular vector factors, go to 50.
344 *
345  IF( icompq.EQ.1 ) THEN
346  go to 50
347  END IF
348 *
349 * The nodes on the bottom level of the tree were solved
350 * by DLASDQ. The corresponding left and right singular vector
351 * matrices are in explicit form. First apply back the left
352 * singular vector matrices.
353 *
354  ndb1 = ( nd+1 ) / 2
355  DO 10 i = ndb1, nd
356 *
357 * IC : center row of each node
358 * NL : number of rows of left subproblem
359 * NR : number of rows of right subproblem
360 * NLF: starting row of the left subproblem
361 * NRF: starting row of the right subproblem
362 *
363  i1 = i - 1
364  ic = iwork( inode+i1 )
365  nl = iwork( ndiml+i1 )
366  nr = iwork( ndimr+i1 )
367  nlf = ic - nl
368  nrf = ic + 1
369  CALL dgemm( 'T', 'N', nl, nrhs, nl, one, u( nlf, 1 ), ldu,
370  $ b( nlf, 1 ), ldb, zero, bx( nlf, 1 ), ldbx )
371  CALL dgemm( 'T', 'N', nr, nrhs, nr, one, u( nrf, 1 ), ldu,
372  $ b( nrf, 1 ), ldb, zero, bx( nrf, 1 ), ldbx )
373  10 continue
374 *
375 * Next copy the rows of B that correspond to unchanged rows
376 * in the bidiagonal matrix to BX.
377 *
378  DO 20 i = 1, nd
379  ic = iwork( inode+i-1 )
380  CALL dcopy( nrhs, b( ic, 1 ), ldb, bx( ic, 1 ), ldbx )
381  20 continue
382 *
383 * Finally go through the left singular vector matrices of all
384 * the other subproblems bottom-up on the tree.
385 *
386  j = 2**nlvl
387  sqre = 0
388 *
389  DO 40 lvl = nlvl, 1, -1
390  lvl2 = 2*lvl - 1
391 *
392 * find the first node LF and last node LL on
393 * the current level LVL
394 *
395  IF( lvl.EQ.1 ) THEN
396  lf = 1
397  ll = 1
398  ELSE
399  lf = 2**( lvl-1 )
400  ll = 2*lf - 1
401  END IF
402  DO 30 i = lf, ll
403  im1 = i - 1
404  ic = iwork( inode+im1 )
405  nl = iwork( ndiml+im1 )
406  nr = iwork( ndimr+im1 )
407  nlf = ic - nl
408  nrf = ic + 1
409  j = j - 1
410  CALL dlals0( icompq, nl, nr, sqre, nrhs, bx( nlf, 1 ), ldbx,
411  $ b( nlf, 1 ), ldb, perm( nlf, lvl ),
412  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
413  $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
414  $ difl( nlf, lvl ), difr( nlf, lvl2 ),
415  $ z( nlf, lvl ), k( j ), c( j ), s( j ), work,
416  $ info )
417  30 continue
418  40 continue
419  go to 90
420 *
421 * ICOMPQ = 1: applying back the right singular vector factors.
422 *
423  50 continue
424 *
425 * First now go through the right singular vector matrices of all
426 * the tree nodes top-down.
427 *
428  j = 0
429  DO 70 lvl = 1, nlvl
430  lvl2 = 2*lvl - 1
431 *
432 * Find the first node LF and last node LL on
433 * the current level LVL.
434 *
435  IF( lvl.EQ.1 ) THEN
436  lf = 1
437  ll = 1
438  ELSE
439  lf = 2**( lvl-1 )
440  ll = 2*lf - 1
441  END IF
442  DO 60 i = ll, lf, -1
443  im1 = i - 1
444  ic = iwork( inode+im1 )
445  nl = iwork( ndiml+im1 )
446  nr = iwork( ndimr+im1 )
447  nlf = ic - nl
448  nrf = ic + 1
449  IF( i.EQ.ll ) THEN
450  sqre = 0
451  ELSE
452  sqre = 1
453  END IF
454  j = j + 1
455  CALL dlals0( icompq, nl, nr, sqre, nrhs, b( nlf, 1 ), ldb,
456  $ bx( nlf, 1 ), ldbx, perm( nlf, lvl ),
457  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
458  $ givnum( nlf, lvl2 ), ldu, poles( nlf, lvl2 ),
459  $ difl( nlf, lvl ), difr( nlf, lvl2 ),
460  $ z( nlf, lvl ), k( j ), c( j ), s( j ), work,
461  $ info )
462  60 continue
463  70 continue
464 *
465 * The nodes on the bottom level of the tree were solved
466 * by DLASDQ. The corresponding right singular vector
467 * matrices are in explicit form. Apply them back.
468 *
469  ndb1 = ( nd+1 ) / 2
470  DO 80 i = ndb1, nd
471  i1 = i - 1
472  ic = iwork( inode+i1 )
473  nl = iwork( ndiml+i1 )
474  nr = iwork( ndimr+i1 )
475  nlp1 = nl + 1
476  IF( i.EQ.nd ) THEN
477  nrp1 = nr
478  ELSE
479  nrp1 = nr + 1
480  END IF
481  nlf = ic - nl
482  nrf = ic + 1
483  CALL dgemm( 'T', 'N', nlp1, nrhs, nlp1, one, vt( nlf, 1 ), ldu,
484  $ b( nlf, 1 ), ldb, zero, bx( nlf, 1 ), ldbx )
485  CALL dgemm( 'T', 'N', nrp1, nrhs, nrp1, one, vt( nrf, 1 ), ldu,
486  $ b( nrf, 1 ), ldb, zero, bx( nrf, 1 ), ldbx )
487  80 continue
488 *
489  90 continue
490 *
491  return
492 *
493 * End of DLALSA
494 *
495  END