LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
slasda.f
Go to the documentation of this file.
1 *> \brief \b SLASDA 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 SLASDA + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasda.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasda.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasda.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLASDA( 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 * REAL 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, SLASDA 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, SLASD0, 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 REAL 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 REAL 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 REAL 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 REAL 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, dimension ( N )
131 *> 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 REAL array, dimension ( LDU, NLVL ),
139 *> where NLVL = floor(log_2 (N/SMLSIZ))).
140 *> \endverbatim
141 *>
142 *> \param[out] DIFR
143 *> \verbatim
144 *> DIFR is REAL 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 SLASD8 for details.
152 *> \endverbatim
153 *>
154 *> \param[out] Z
155 *> \verbatim
156 *> Z is REAL 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 REAL 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, dimension ( LDGCOL, NLVL )
201 *> 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 REAL 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 REAL 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 REAL 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 REAL 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 *> \date September 2012
262 *
263 *> \ingroup auxOTHERauxiliary
264 *
265 *> \par Contributors:
266 * ==================
267 *>
268 *> Ming Gu and Huan Ren, Computer Science Division, University of
269 *> California at Berkeley, USA
270 *>
271 * =====================================================================
272  SUBROUTINE slasda( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K,
273  $ difl, difr, z, poles, givptr, givcol, ldgcol,
274  $ perm, givnum, c, s, work, iwork, info )
275 *
276 * -- LAPACK auxiliary routine (version 3.4.2) --
277 * -- LAPACK is a software package provided by Univ. of Tennessee, --
278 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
279 * September 2012
280 *
281 * .. Scalar Arguments ..
282  INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE
283 * ..
284 * .. Array Arguments ..
285  INTEGER GIVCOL( ldgcol, * ), GIVPTR( * ), IWORK( * ),
286  $ k( * ), perm( ldgcol, * )
287  REAL C( * ), D( * ), DIFL( ldu, * ), DIFR( ldu, * ),
288  $ e( * ), givnum( ldu, * ), poles( ldu, * ),
289  $ s( * ), u( ldu, * ), vt( ldu, * ), work( * ),
290  $ z( ldu, * )
291 * ..
292 *
293 * =====================================================================
294 *
295 * .. Parameters ..
296  REAL ZERO, ONE
297  parameter( zero = 0.0e+0, one = 1.0e+0 )
298 * ..
299 * .. Local Scalars ..
300  INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK,
301  $ j, lf, ll, lvl, lvl2, m, ncc, nd, ndb1, ndiml,
302  $ ndimr, nl, nlf, nlp1, nlvl, nr, nrf, nrp1, nru,
303  $ nwork1, nwork2, smlszp, sqrei, vf, vfi, vl, vli
304  REAL ALPHA, BETA
305 * ..
306 * .. External Subroutines ..
307  EXTERNAL scopy, slasd6, slasdq, slasdt, slaset, xerbla
308 * ..
309 * .. Executable Statements ..
310 *
311 * Test the input parameters.
312 *
313  info = 0
314 *
315  IF( ( icompq.LT.0 ) .OR. ( icompq.GT.1 ) ) THEN
316  info = -1
317  ELSE IF( smlsiz.LT.3 ) THEN
318  info = -2
319  ELSE IF( n.LT.0 ) THEN
320  info = -3
321  ELSE IF( ( sqre.LT.0 ) .OR. ( sqre.GT.1 ) ) THEN
322  info = -4
323  ELSE IF( ldu.LT.( n+sqre ) ) THEN
324  info = -8
325  ELSE IF( ldgcol.LT.n ) THEN
326  info = -17
327  END IF
328  IF( info.NE.0 ) THEN
329  CALL xerbla( 'SLASDA', -info )
330  RETURN
331  END IF
332 *
333  m = n + sqre
334 *
335 * If the input matrix is too small, call SLASDQ to find the SVD.
336 *
337  IF( n.LE.smlsiz ) THEN
338  IF( icompq.EQ.0 ) THEN
339  CALL slasdq( 'U', sqre, n, 0, 0, 0, d, e, vt, ldu, u, ldu,
340  $ u, ldu, work, info )
341  ELSE
342  CALL slasdq( 'U', sqre, n, m, n, 0, d, e, vt, ldu, u, ldu,
343  $ u, ldu, work, info )
344  END IF
345  RETURN
346  END IF
347 *
348 * Book-keeping and set up the computation tree.
349 *
350  inode = 1
351  ndiml = inode + n
352  ndimr = ndiml + n
353  idxq = ndimr + n
354  iwk = idxq + n
355 *
356  ncc = 0
357  nru = 0
358 *
359  smlszp = smlsiz + 1
360  vf = 1
361  vl = vf + m
362  nwork1 = vl + m
363  nwork2 = nwork1 + smlszp*smlszp
364 *
365  CALL slasdt( n, nlvl, nd, iwork( inode ), iwork( ndiml ),
366  $ iwork( ndimr ), smlsiz )
367 *
368 * for the nodes on bottom level of the tree, solve
369 * their subproblems by SLASDQ.
370 *
371  ndb1 = ( nd+1 ) / 2
372  DO 30 i = ndb1, nd
373 *
374 * IC : center row of each node
375 * NL : number of rows of left subproblem
376 * NR : number of rows of right subproblem
377 * NLF: starting row of the left subproblem
378 * NRF: starting row of the right subproblem
379 *
380  i1 = i - 1
381  ic = iwork( inode+i1 )
382  nl = iwork( ndiml+i1 )
383  nlp1 = nl + 1
384  nr = iwork( ndimr+i1 )
385  nlf = ic - nl
386  nrf = ic + 1
387  idxqi = idxq + nlf - 2
388  vfi = vf + nlf - 1
389  vli = vl + nlf - 1
390  sqrei = 1
391  IF( icompq.EQ.0 ) THEN
392  CALL slaset( 'A', nlp1, nlp1, zero, one, work( nwork1 ),
393  $ smlszp )
394  CALL slasdq( 'U', sqrei, nl, nlp1, nru, ncc, d( nlf ),
395  $ e( nlf ), work( nwork1 ), smlszp,
396  $ work( nwork2 ), nl, work( nwork2 ), nl,
397  $ work( nwork2 ), info )
398  itemp = nwork1 + nl*smlszp
399  CALL scopy( nlp1, work( nwork1 ), 1, work( vfi ), 1 )
400  CALL scopy( nlp1, work( itemp ), 1, work( vli ), 1 )
401  ELSE
402  CALL slaset( 'A', nl, nl, zero, one, u( nlf, 1 ), ldu )
403  CALL slaset( 'A', nlp1, nlp1, zero, one, vt( nlf, 1 ), ldu )
404  CALL slasdq( 'U', sqrei, nl, nlp1, nl, ncc, d( nlf ),
405  $ e( nlf ), vt( nlf, 1 ), ldu, u( nlf, 1 ), ldu,
406  $ u( nlf, 1 ), ldu, work( nwork1 ), info )
407  CALL scopy( nlp1, vt( nlf, 1 ), 1, work( vfi ), 1 )
408  CALL scopy( nlp1, vt( nlf, nlp1 ), 1, work( vli ), 1 )
409  END IF
410  IF( info.NE.0 ) THEN
411  RETURN
412  END IF
413  DO 10 j = 1, nl
414  iwork( idxqi+j ) = j
415  10 CONTINUE
416  IF( ( i.EQ.nd ) .AND. ( sqre.EQ.0 ) ) THEN
417  sqrei = 0
418  ELSE
419  sqrei = 1
420  END IF
421  idxqi = idxqi + nlp1
422  vfi = vfi + nlp1
423  vli = vli + nlp1
424  nrp1 = nr + sqrei
425  IF( icompq.EQ.0 ) THEN
426  CALL slaset( 'A', nrp1, nrp1, zero, one, work( nwork1 ),
427  $ smlszp )
428  CALL slasdq( 'U', sqrei, nr, nrp1, nru, ncc, d( nrf ),
429  $ e( nrf ), work( nwork1 ), smlszp,
430  $ work( nwork2 ), nr, work( nwork2 ), nr,
431  $ work( nwork2 ), info )
432  itemp = nwork1 + ( nrp1-1 )*smlszp
433  CALL scopy( nrp1, work( nwork1 ), 1, work( vfi ), 1 )
434  CALL scopy( nrp1, work( itemp ), 1, work( vli ), 1 )
435  ELSE
436  CALL slaset( 'A', nr, nr, zero, one, u( nrf, 1 ), ldu )
437  CALL slaset( 'A', nrp1, nrp1, zero, one, vt( nrf, 1 ), ldu )
438  CALL slasdq( 'U', sqrei, nr, nrp1, nr, ncc, d( nrf ),
439  $ e( nrf ), vt( nrf, 1 ), ldu, u( nrf, 1 ), ldu,
440  $ u( nrf, 1 ), ldu, work( nwork1 ), info )
441  CALL scopy( nrp1, vt( nrf, 1 ), 1, work( vfi ), 1 )
442  CALL scopy( nrp1, vt( nrf, nrp1 ), 1, work( vli ), 1 )
443  END IF
444  IF( info.NE.0 ) THEN
445  RETURN
446  END IF
447  DO 20 j = 1, nr
448  iwork( idxqi+j ) = j
449  20 CONTINUE
450  30 CONTINUE
451 *
452 * Now conquer each subproblem bottom-up.
453 *
454  j = 2**nlvl
455  DO 50 lvl = nlvl, 1, -1
456  lvl2 = lvl*2 - 1
457 *
458 * Find the first node LF and last node LL on
459 * the current level LVL.
460 *
461  IF( lvl.EQ.1 ) THEN
462  lf = 1
463  ll = 1
464  ELSE
465  lf = 2**( lvl-1 )
466  ll = 2*lf - 1
467  END IF
468  DO 40 i = lf, ll
469  im1 = i - 1
470  ic = iwork( inode+im1 )
471  nl = iwork( ndiml+im1 )
472  nr = iwork( ndimr+im1 )
473  nlf = ic - nl
474  nrf = ic + 1
475  IF( i.EQ.ll ) THEN
476  sqrei = sqre
477  ELSE
478  sqrei = 1
479  END IF
480  vfi = vf + nlf - 1
481  vli = vl + nlf - 1
482  idxqi = idxq + nlf - 1
483  alpha = d( ic )
484  beta = e( ic )
485  IF( icompq.EQ.0 ) THEN
486  CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
487  $ work( vfi ), work( vli ), alpha, beta,
488  $ iwork( idxqi ), perm, givptr( 1 ), givcol,
489  $ ldgcol, givnum, ldu, poles, difl, difr, z,
490  $ k( 1 ), c( 1 ), s( 1 ), work( nwork1 ),
491  $ iwork( iwk ), info )
492  ELSE
493  j = j - 1
494  CALL slasd6( icompq, nl, nr, sqrei, d( nlf ),
495  $ work( vfi ), work( vli ), alpha, beta,
496  $ iwork( idxqi ), perm( nlf, lvl ),
497  $ givptr( j ), givcol( nlf, lvl2 ), ldgcol,
498  $ givnum( nlf, lvl2 ), ldu,
499  $ poles( nlf, lvl2 ), difl( nlf, lvl ),
500  $ difr( nlf, lvl2 ), z( nlf, lvl ), k( j ),
501  $ c( j ), s( j ), work( nwork1 ),
502  $ iwork( iwk ), info )
503  END IF
504  IF( info.NE.0 ) THEN
505  RETURN
506  END IF
507  40 CONTINUE
508  50 CONTINUE
509 *
510  RETURN
511 *
512 * End of SLASDA
513 *
514  END
subroutine slasd6(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)
SLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by...
Definition: slasd6.f:316
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slasdt(N, LVL, ND, INODE, NDIML, NDIMR, MSUB)
SLASDT creates a tree of subproblems for bidiagonal divide and conquer. Used by sbdsdc.
Definition: slasdt.f:107
subroutine slasdq(UPLO, SQRE, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e...
Definition: slasdq.f:213
subroutine slasda(ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, PERM, GIVNUM, C, S, WORK, IWORK, INFO)
SLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
Definition: slasda.f:275
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112