LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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, dimension (7*N)
243 *> \endverbatim
244 *>
245 *> \param[out] INFO
246 *> \verbatim
247 *> INFO is INTEGER
248 *> = 0: successful exit.
249 *> < 0: if INFO = -i, the i-th argument had an illegal value.
250 *> > 0: if INFO = 1, a singular value did not converge
251 *> \endverbatim
252 *
253 * Authors:
254 * ========
255 *
256 *> \author Univ. of Tennessee
257 *> \author Univ. of California Berkeley
258 *> \author Univ. of Colorado Denver
259 *> \author NAG Ltd.
260 *
261 *> \ingroup OTHERauxiliary
262 *
263 *> \par Contributors:
264 * ==================
265 *>
266 *> Ming Gu and Huan Ren, Computer Science Division, University of
267 *> California at Berkeley, USA
268 *>
269 * =====================================================================
270  SUBROUTINE dlasda( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
271  $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL,
272  $ PERM, GIVNUM, C, S, WORK, IWORK, INFO )
273 *
274 * -- LAPACK auxiliary routine --
275 * -- LAPACK is a software package provided by Univ. of Tennessee, --
276 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
277 *
278 * .. Scalar Arguments ..
279  INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
280 * ..
281 * .. Array Arguments ..
282  INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ),
283  $ K( * ), PERM( LDGCOL, * )
284  DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ),
285  $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ),
286  $ s( * ), u( ldu, * ), vt( ldu, * ), work( * ),
287  $ z( ldu, * )
288 * ..
289 *
290 * =====================================================================
291 *
292 * .. Parameters ..
293  DOUBLE PRECISION ZERO, ONE
294  PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
295 * ..
296 * .. Local Scalars ..
297  INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
298  $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML,
299  $ ndimr, nl, nlf, nlp1, nlvl, nr, nrf, nrp1, nru,
300  $ nwork1, nwork2, smlszp, sqrei, vf, vfi, vl, vli
301  DOUBLE PRECISION ALPHA, BETA
302 * ..
303 * .. External Subroutines ..
304  EXTERNAL dcopy, dlasd6, dlasdq, dlasdt, dlaset, xerbla
305 * ..
306 * .. Executable Statements ..
307 *
308 * Test the input parameters.
309 *
310  info = 0
311 *
312  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
313  info = -1
314  ELSE IF( smlsiz.LT.3 ) THEN
315  info = -2
316  ELSE IF( n.LT.0 ) THEN
317  info = -3
318  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
319  info = -4
320  ELSE IF( ldu.LT.( n+sqre ) ) THEN
321  info = -8
322  ELSE IF( ldgcol.LT.n ) THEN
323  info = -17
324  END IF
325  IF( info.NE.0 ) THEN
326  CALL xerbla( 'DLASDA', -info )
327  RETURN
328  END IF
329 *
330  m = n + sqre
331 *
332 * If the input matrix is too small, call DLASDQ to find the SVD.
333 *
334  IF( n.LE.smlsiz ) THEN
335  IF( icompq.EQ.0 ) THEN
336  CALL dlasdq( 'U', sqre, n, 0, 0, 0, d, e, vt, ldu, u, ldu,
337  $ u, ldu, work, info )
338  ELSE
339  CALL dlasdq( 'U', sqre, n, m, n, 0, d, e, vt, ldu, u, ldu,
340  $ u, ldu, work, info )
341  END IF
342  RETURN
343  END IF
344 *
345 * Book-keeping and set up the computation tree.
346 *
347  inode = 1
348  ndiml = inode + n
349  ndimr = ndiml + n
350  idxq = ndimr + n
351  iwk = idxq + n
352 *
353  ncc = 0
354  nru = 0
355 *
356  smlszp = smlsiz + 1
357  vf = 1
358  vl = vf + m
359  nwork1 = vl + m
360  nwork2 = nwork1 + smlszp*smlszp
361 *
362  CALL dlasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
363  $ iwork( ndimr ), smlsiz )
364 *
365 * for the nodes on bottom level of the tree, solve
366 * their subproblems by DLASDQ.
367 *
368  ndb1 = ( nd+1 ) / 2
369  DO 30 i = ndb1, nd
370 *
371 * IC : center row of each node
372 * NL : number of rows of left subproblem
373 * NR : number of rows of right subproblem
374 * NLF: starting row of the left subproblem
375 * NRF: starting row of the right subproblem
376 *
377  i1 = i - 1
378  ic = iwork( inode+i1 )
379  nl = iwork( ndiml+i1 )
380  nlp1 = nl + 1
381  nr = iwork( ndimr+i1 )
382  nlf = ic - nl
383  nrf = ic + 1
384  idxqi = idxq + nlf - 2
385  vfi = vf + nlf - 1
386  vli = vl + nlf - 1
387  sqrei = 1
388  IF( icompq.EQ.0 ) THEN
389  CALL dlaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
390  $ smlszp )
391  CALL dlasdq( 'U', sqrei, nl, nlp1, nru, ncc, d( nlf ),
392  $ e( nlf ), work( nwork1 ), smlszp,
393  $ work( nwork2 ), nl, work( nwork2 ), nl,
394  $ work( nwork2 ), info )
395  itemp = nwork1 + nl*smlszp
396  CALL dcopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
397  CALL dcopy( nlp1, work( itemp ), 1, work( vli ), 1 )
398  ELSE
399  CALL dlaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
400  CALL dlaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ), ldu )
401  CALL dlasdq( 'U', sqrei, nl, nlp1, nl, ncc, d( nlf ),
402  $ e( nlf ), vt( nlf, 1 ), ldu, u( nlf, 1 ), ldu,
403  $ u( nlf, 1 ), ldu, work( nwork1 ), info )
404  CALL dcopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
405  CALL dcopy( nlp1, vt( nlf, nlp1 ), 1, work( vli ), 1 )
406  END IF
407  IF( info.NE.0 ) THEN
408  RETURN
409  END IF
410  DO 10 j = 1, nl
411  iwork( idxqi+j ) = j
412  10 CONTINUE
413  IF( ( i.EQ.nd ) .AND. ( sqre.EQ.0 ) ) THEN
414  sqrei = 0
415  ELSE
416  sqrei = 1
417  END IF
418  idxqi = idxqi + nlp1
419  vfi = vfi + nlp1
420  vli = vli + nlp1
421  nrp1 = nr + sqrei
422  IF( icompq.EQ.0 ) THEN
423  CALL dlaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
424  $ smlszp )
425  CALL dlasdq( 'U', sqrei, nr, nrp1, nru, ncc, d( nrf ),
426  $ e( nrf ), work( nwork1 ), smlszp,
427  $ work( nwork2 ), nr, work( nwork2 ), nr,
428  $ work( nwork2 ), info )
429  itemp = nwork1 + ( nrp1-1 )*smlszp
430  CALL dcopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
431  CALL dcopy( nrp1, work( itemp ), 1, work( vli ), 1 )
432  ELSE
433  CALL dlaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
434  CALL dlaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ), ldu )
435  CALL dlasdq( 'U', sqrei, nr, nrp1, nr, ncc, d( nrf ),
436  $ e( nrf ), vt( nrf, 1 ), ldu, u( nrf, 1 ), ldu,
437  $ u( nrf, 1 ), ldu, work( nwork1 ), info )
438  CALL dcopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
439  CALL dcopy( nrp1, vt( nrf, nrp1 ), 1, work( vli ), 1 )
440  END IF
441  IF( info.NE.0 ) THEN
442  RETURN
443  END IF
444  DO 20 j = 1, nr
445  iwork( idxqi+j ) = j
446  20 CONTINUE
447  30 CONTINUE
448 *
449 * Now conquer each subproblem bottom-up.
450 *
451  j = 2**nlvl
452  DO 50 lvl = nlvl, 1, -1
453  lvl2 = lvl*2 - 1
454 *
455 * Find the first node LF and last node LL on
456 * the current level LVL.
457 *
458  IF( lvl.EQ.1 ) THEN
459  lf = 1
460  ll = 1
461  ELSE
462  lf = 2**( lvl-1 )
463  ll = 2*lf - 1
464  END IF
465  DO 40 i = lf, ll
466  im1 = i - 1
467  ic = iwork( inode+im1 )
468  nl = iwork( ndiml+im1 )
469  nr = iwork( ndimr+im1 )
470  nlf = ic - nl
471  nrf = ic + 1
472  IF( i.EQ.ll ) THEN
473  sqrei = sqre
474  ELSE
475  sqrei = 1
476  END IF
477  vfi = vf + nlf - 1
478  vli = vl + nlf - 1
479  idxqi = idxq + nlf - 1
480  alpha = d( ic )
481  beta = e( ic )
482  IF( icompq.EQ.0 ) THEN
483  CALL dlasd6( icompq, nl, nr, sqrei, d( nlf ),
484  $ work( vfi ), work( vli ), alpha, beta,
485  $ iwork( idxqi ), perm, givptr( 1 ), givcol,
486  $ ldgcol, givnum, ldu, poles, difl, difr, z,
487  $ k( 1 ), c( 1 ), s( 1 ), work( nwork1 ),
488  $ iwork( iwk ), info )
489  ELSE
490  j = j - 1
491  CALL dlasd6( icompq, nl, nr, sqrei, d( nlf ),
492  $ work( vfi ), work( vli ), alpha, beta,
493  $ iwork( idxqi ), perm( nlf, lvl ),
494  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
495  $ givnum( nlf, lvl2 ), ldu,
496  $ poles( nlf, lvl2 ), difl( nlf, lvl ),
497  $ difr( nlf, lvl2 ), z( nlf, lvl ), k( j ),
498  $ c( j ), s( j ), work( nwork1 ),
499  $ iwork( iwk ), info )
500  END IF
501  IF( info.NE.0 ) THEN
502  RETURN
503  END IF
504  40 CONTINUE
505  50 CONTINUE
506 *
507  RETURN
508 *
509 * End of DLASDA
510 *
511  END
subroutine dlasd6(ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, IWORK, INFO)
DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by...
Definition: dlasd6.f:313
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine dlasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
DLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
Definition: dlasdt.f:105
subroutine dlasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: dlasda.f:273
subroutine dlasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
Definition: dlasdq.f:211
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82