LAPACK  3.6.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.
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
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:112
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:213
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:276
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:107
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:316
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53