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