LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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